AD-A211  918 


MASSACHUSETTS  INSTITUTE  OF  TECHNOLOGY 


VLSI  PUBLICATIONS 


VLSI  Memo  No.  89-529 
May  1989 


nn 


i J  J  £ 


I.  „ 


It  F  ECT  £  SSfc 
SEP  0  5 1989 1  | 

0 b 


Sequential  Design  of  Experiments  with  Physically  Based  Models 

Michele  E.  Storm 


Abstract 


The  application  of  physically  based  models  to  sequential  optimization  was  explored  and  the 
benefits  measured  by  comparison  to  optimizations  performed  using  control  variable 
polynomial  models.  The  optimization  algorithm  developed  is  based  on  the  sequential  use 
of  "local"  (weighted)  linear  regression  models.  A  new  operating  point  or  design  is 
recommended  at  the  optimum  of  the  model  within  the  region  where  the  model  is 
considered  valid.  This  region,  within  which  extrapolations  based  on  the  model  are  believed 
to  be  sufficiently  accurate,  is  defined  by  constraints  based  on  the  estimated  predictive  error 
of  the  model  and  the  distance  from  the  data.  A  new  model  is  created  after  each  data  point 
is  collected. 


As  a  test  case,  the  sequential  optimizer  was  applied  to  the  design  of  a  paper  helicopter  for 
maximum  time  of  flight.  The  physically  based  model  of  the  paper  helicopter  was 
developed  through  the  use  of  dimensional  analysis,  a  technique  which  groups  variables 
according  to  their  dimensions.  It  was  found  that  the  physically  based  model  improved  the 
design  of  the  helicopter  more  rapidly  than  the  polynomial  models.  For  example,  in  one 
comparison  of  the  physical  model  and  the  linear  model,  the  physical  model  reached  the 
optimum  design  after  four  sequential  designs,  while  the  linear  model  hadn’t  reached  the 
optimum  after  ten  designs.  . ,  _ 


The  superior  performance  of  the  optimizer  with  the  physically  based  model  is  attributed  to 
the  model  providing  a  truer  representation  of  the  actual  air-helicopter  system,  resulting  in  a 
larger  region  where  the  model  is  valid  and  a  better  choice  of  direction  within  that  region. 
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The  application  of  physically  based  models  to  sequential  optimization  was  explored  and  the 
benefits  measured  by  comparison  to  optimizations  performed  using  control  variable  polynomial 
models.  The  optimization  algorithm  developed  is  based  on  the  sequential  use  of  "local"  (weighted) 
linear  regression  models.  A  new  operating  point  or  design  is  recommended  at  the  optimum  of  the 
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based  on  the  estimated  predictive  error  of  the  model  and  the  distance  from  the  data.  A  new  model 
is  created  after  each  data  point  is  collected. 

As  a  test  case,  the  sequential  optimizer  was  applied  to  the  design  of  a  paper  helicopter  for 
maximum  time  of  flight.  The  physically  based  model  of  the  paper  helicopter  was  developed 
through  the  use  of  dimensional  analysis,  a  technique  which  groups  variables  according  to  their 
dimensions.  It  was  found  that  the  physically  based  model  improved  the  design  of  the  helicopter 
more  rapidly  than  the  polynomial  models.  For  example,  in  one  comparison  of  the  physical  model 
and  the  linear  model,  the  physical  model  reached  the  optimum  design  after  four  sequential  designs, 
while  the  linear  model  hadn’t  reached  the  optimum  after  ten  designs. 

The  superior  performance  of  the  optimizer  with  the  physically  based  model  is  attributed  to 
the  model  providing  a  truer  representation  of  the  actual  air-helicopter  system,  resulting  in  a  larger 
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List  of  Svmbols 


y:  the  objective  function  to  be  maximized  by  the  sequential  optimizer 

x:  the  vector  of  control  variables  -  the  actual  variables  which  the  experimenter  changes 

z:  the  vector  of  predictor  variables  in  the  regression  function.  The  elements  of  z  are  a 
function  of  the  elements  of  x. 

Z:  the  design  matrix  whose  rows  are  z* 

C,:  vector  of  random  error  terms  associated  with  the  linear  model 

V:  the  matrix  which  contains  the  terms,  w,  that  determine  the  relative  magnitude  of  the 
error  variance  with  respect  to  the  standard  error  variance 

o2:  the  standard  error  variance,  the  variance  of  the  reference  point. 

w:  the  the  weighting  factor  for  the  error  variance 

A;:  scaling  constants  which  determine  the  local  region  by  determining  w 

Q:  matrix  which  contains  A,  and  is  part  of  the  weighting  function 

P  :  the  vector  of  coefficients  in  the  linear  regression  model 

s2:  the  estimated  standard  variance 

Z:  the  covariance  matrix  of  x 

p:  the  mean  vector  of  x 

L:  the  length  of  the  paper  helicopter  blade 

b:  the  width  of  the  paper  helicopter  wing 

wt:  the  total  weight  of  the  paper  helicopter 

wa :  the  weight  of  the  paper  added  to  the  end  of  the  helicopter  tail 

wp:  the  weight  of  the  paper  in  the  helicopter 

V:  the  descent  velocity  of  the  helicopter 
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1  INTRODUCTION 


1.1  Motivation 

One  of  the  fundamental  goals  of  methods  of  experimental  design  is  the  realization  of 
improvement  of  products  and  processes  with  a  minimum  of  experimental  effort.  The 
work  described  in  this  thesis  seeks  to  minimize  the  number  of  experiments  needed  to 
optimize  the  design  of  a  product  or  operation  of  a  process  through  the  use  of  physically- 
based  models  in  an  algorithm  based  on  sequential  design  of  experiments. 

1.2  Related  Work 

There  are  two  primary  methods  of  experimental  design;  parallel  and  sequential. 
Both  methods  are  well  known  and  have  been  written  about  extensively  in  the  literature. 
In  parallel  design  of  experiments  (DOE),  a  set  of  experiments  is  run  and  then  the  results 
are  used  to  optimize  the  product.  Classical  DOE,  such  as  partial  and  full  factorial  designs 
(Box,  Hunter,  and  Hunter,  1978)  and  Box-Behnken  designs  (Box  and  Behnken,1960), 
as  well  as  the  orthogonal  arrays  popular  in  quality  control  (Taguchi,  1986)  fall  into  this 
category . 

In  sequential  DOE  data  are  analyzed  as  they  are  collected,  and  used  to  determine  the 
next  experiment  to  be  run.  The  advantage  of  this  scheme  is  that  it  allows  the  the 
experimenter  to  use  the  information  available  in  the  collected  data  to  tailor  the  remaining 
experiments  to  more  effectively  optimize  the  process.  Examples  of  this  method  of 
experimentation  are  the  Simplex  method  and  the  Ultramax  method  (Moreno,  1986). 

The  key  issue  in  deciding  if  sequentially  designed  experiments  are  applicable  is  the 
speed  with  which  feedback  from  the  experiments  can  be  obtained.  For  example,  in 
agricultural  experiments,  which  require  an  entire  growing  season  to  complete,  parallel 
experiments  make  are  the  most  suitable  because  a  series  of  sequentially  performed 
experiments  would  require  a  prohibitively  long  time  to  complete.  On  the  other  hand, 
when  experiments  are  performed  on  manufacturing  production  lines  the  product  is 
usually  produced  sequentially  and  measurements  on  the  response  of  interest,  such  as 
the  dimensions  of  a  part,  can  be  obtained  almost  immediately.  In  this  case, 
sequentially  designed  experiments  are  most  appropriate. 
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The  concept  of  on-line  process  optimization  through  repeated  use  of  designed 
experiments  was  first  introduced  by  G.  E.  P.  Box  in  "Evolutionary  Operation”.  This 
idea  was  further  developed  by  the  Ultramax  Corporation,  which  currently  markets  a 
computer  program  to  perform  sequential  opdmization.  This  software  performs  sequential 
optimization  based  on  quadratic  polynomials  relating  the  response  to  the  predictor 
variables.  Our  work  builds  on  these  sources. 

1.3  Physical  Models 

In  many  engineering  applications  of  statistical  methods,  including  the  Ultramax 
Method,  a  physical  system  is  modeled  by  a  polynomial.  This  approach  assumes  that  a 
polynomial  can  adequately  describe  the  relationship  between  the  response  and  the 
predictor  variables.  However,  the  engineer  often  has  some  idea  about  the  relationship 
between  the  response  and  the  predictor  variables,  for  example,  in  the  case  of  thin  film 
deposition  process,  it  may  be  known  that  the  growth  rate  varies  inversely  with  the 
pressure  (P)  and  proportionally  to  flow  rate  (Q).  In  this  case  a  form  of  the  model  that 
relates  the  grouped  variable  Q/P  to  growth  rate  would  be  more  accurate  than  one  that 
related  P  and  Q  to  growth  rate.  Using  the  engineers'  prior  knowledge  of  a  process  results 
in  the  choice  of  a  model  form  which  more  accurately  represents  the  true  situation.  Both 
physically  based  modelling  and  sequential  DOE  deal  with  the  efficient  incorporation  of  the 
available  information. 

1.4  Current  Work 

Our  immediate  objective  is  to  determine  how  physically  based  models  enhance  the 
performance  of  sequential  optimizers.  Our  long-term  objective  is  to  explore  the 
optimization  of  manufacturing  processes  using  dedicated  sequential  optimization 
algorithms  with  embedded  physically  based  models.  The  type  of  physically  based  model 
examined  in  this  thesis  is  derived  from  what  is  known  in  the  physical  sciences  as 
dimensional  analysis. 

In  order  to  achieve  our  goal  a  very  basic  sequential  optimizer  was  developed.  This 
optimizer  drew  on  the  concepts  of  the  Ultramax  Method.  Our  code  has  the  flexibility  to 
allow  the  use  of  any  model  which  is  linear  in  the  estimated  coefficients.  This  flexibility 
was  needed  to  implement  some  of  the  models.  Also,  developing  our  own  simple 
sequential  optimizer  has  the  advantage  that  we  know  the  algorithms  involved  so  tha  the 


6 


effect  of  the  physically  based  models  can  be  more  precisely  evaluated.  The  design  of  a 
paper  helicopter  was  used  to  test  the  ability  of  physical  models  to  improve  the  performance 
of  a  sequential  optimizer. 
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2  DE\  LLOPMENT  OF  A  SEQUENTIAL  OPTIMIZER 


The  sequential  optimizer  developed  for  this  thesis  works  in  the  following  manner. 
Based  on  the  available  data  the  optimizer  suggests  a  set  of  control  variables  which  it 
believes  to  be  an  improvement  over  the  last  run.  The  process  operator  then  runs  the 
process  at  these  settings.  The  quantity  of  interest  is  recorded  as  well  as  the  settings  at 
which  the  process  was  run.  This  data  is  then  added  to  the  data  file  which  the  optimizer 
uses  to  determine  the  next  point  that  it  will  suggest.  The  cycle  of  running  the  code  to  obtain 
a  suggested  new'  set  of  control  variables  at  which  to  run  then  repeats  itself  (see  Figure  1). 
In  this  manner  the  process  is  gradually  improved. 


Figure  1:  Diagram  of  process  with  sequential  optimizer 

In  the  code  developed  for  this  thesis  there  are  three  steps  involved  in  determining 
the  next  recommended  operating  point.  The  first  step  is  the  creation  of  a  model  relating  the 
input  parameters  to  the  objective  function.  The  second  step  is  the  determination  of  the 
region  in  which  the  model  is  an  adequate  representation  of  the  relationship  between  the 
objective  function  and  the  control  variables.  The  last  step  is  the  optimization  of  the 
objective  function  based  on  the  model  created  in  step  1.  The  solution  is  constrained  to  be  in 
the  region  in  which  the  model  is  considered  valid  as  determined  by  step  2.  The 
optimization  is  also  subject  to  the  constraints  that  the  user  places  on  the  control  variables. 


8 


2.1  Model  and  Notation 


The  first  basic  assumption  is  that  for  each  individual  response  the  relationship  between  the 
objective  function  (y)  and  the  control  variables  (x)  can  be  locally  described  by  an  equation 
of  the  form 


y  =  Z(3  +  C 


(1) 


where 

C,  ~  N(0.Vc2)  and  z  =  f(x) 


(2) 


The  V  term  is  defined  by 


w'here 


[  Wl  O' 

W2 

.  0  wn  _ 


(3) 


w,  =  expUxj  -  xr)tQ(xi  -  xr)/(2m)j 


(4) 


and 


Q  = 


o 


(5) 


The  reference  point  around  which  the  model  is  created  is  denoted  xr  and  m  is  the  number  of 
control  variables.  The  diagonal  elements  of  Q  are  a  function  of  Aj  which  is  the  length  of 
the  local  region  for  the  control  variable  Xj. 

There  are  three  important  features  of  this  formulation.  The  magnitude  of  the  error, 
Q,  and  hence  its  variance,  WjO2  =  exp[(Xj  -  x^K^Cxj  -  xr)/(2m))a2,  is  a  function  of  x  and 
Q.  The  predictor  variables,  z,  are  a  function  of  x,  and  the  model  is  linear  in  the 
coefficients,  P. 
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The  first  important  feature  is  the  treatment  of  the  error  variance  as  a  function  of  x. 
In  the  standard  formulation  of  the  linear  regression  problem,  the  error  variance  is  assumed 
to  be  the  same  for  every  data  point.  (The  V  in  the  present  error  distribution  is  replaced  by 
the  identity  matrix  I  in  the  standard  formulation).  In  our  formulation,  the  variance  of  the 
error  increases  exponentially  as  (x^xrl^Kxi-Xr).  The  quantity  (Xi-x^KjIXi-Xr)  can  be  seen 
to  represent  the  square  of  the  scaled  distance  from  xr  by  performing  the  matrix 
multiplication  of  the  quantity 


(x,-xr)tQ(xi-xr) 


yielding: 


(xirxrl)VAi  +  (xl2-xr2)2/A2  +•••  +  (xim-xrm)2/Ar 


(6) 

(7) 


Q  makes  the  distance  dimensionless  so  that  variables  with  different  units  can  be  compared. 
Most  importantly,  it  defines  the  size  and  shape  of  the  local  region  by  determining  the  rate  at 
which  movement  away  from  the  reference  point  in  a  particular  variable  is  penalized. 

The  use  of  the  exponentially  increasing  error  reflects  the  assumption  that  die  model 
is  insufficient  to  describe  the  process  well  over  the  entire  range  of  data  but  that  it  can  give 
an  adequate  description  of  the  process  in  a  smaller  region.  The  model  is  created  to  be  most 
accurate  near  the  reference  point.  Inadequacies  in  the  model  are  treated  as  if  they  were  pure 
error.  Although  this  is  incorrect,  it  provides  a  useful  way  to  treat  model  error.  The  fact 
that  the  error  variance  isn't  constant  affects  the  way  that  the  coefficients  (5  are  estimated,  as 
well  as  the  error  in  the  prediction  of  y. 


Another  feature  of  this  model  is  that  y  is  a  function  of  z,  which  is  in  turn  a  function 
of  x.  The  vector  z  is  used  to  represent  the  transformation  of  the  data  from  the  original 
control  variables.  For  example,  if,  in  terms  of  the  original  control  variables,  the  model  is 

y  =  P0  +  Pi  Vx7/x3  +  P2x2  (8) 

then 

Zq  —  1,  Zj  =  v  x  j  /X3 ,  Z2  —  X2.  (9) 
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Although  this  is  largely  a  matter  of  notational  convenience,  it  helps  to  clarify  the  fact 
that  the  model  of  the  system  is  more  complex  than  a  simple  linear  relationship  between  y 
and  x.  There  are  many  choices  for  the  relationship  between  z  and  x.  One  of  the  major 
focuses  of  this  thesis  is  the  judicious  choice  of  transformations  from  x  to  z. 

The  third  feature  is  that  the  model  is  linear  in  the  coefficients.  This  means  that  a 
least  squared  error  regression  can  be  used  to  estimate  the  coefficients.  However,  because 
the  error  variance  changes  from  one  data  point  to  another  a  weighted  least  squares 
regression  is  the  appropriate  regression  type  in  this  case.  A  weighted  least  squares 
regression  minimizes: 


Z  wjHy.-y,)2  (10) 

i 

while  an  unweighted  regression  minimizes: 

X(y.-y.)2-  (ID 

i 

The  standard  equation  for  a  weighted  least  squares  regression  is: 

p  =  (ZlV-1Z)'1ZlV1y  (12) 

The  resulting  estimate  of  P  is  distributed  as 

✓v 

P  ~  NCP.CZ'V^Z)*1).  (13) 

Brief  derivations  of  equations  (12)  and  (13)  are  given  in  Appendix  A. 

2.2  Determination  of  the  Region  in  which  the  Model  is  Valid 

In  order  to  determine  the  region  in  which  the  model  will  be  considered  valid,  two 
criteria  were  used.  The  first  is  based  on  the  prediction  error  of  y.  The  second  criteria  is 
based  on  the  "distance"  from  the  data  that  created  the  model.  The  first  criterion  is 
statistically  based.  If  the  assumptions  about  the  form  of  the  model  and  the  nature  of  the 
variance  are  correct,  then  the  90%  confidence  interval  for  the  prediction  value  is  given  by: 

y0  =  yo  ±  t.osV  z^Z^Z^zc^+wos2  (14) 


11 


or 


y0  =  ?0  ±  t.05^  *b(ZlVlZ)  z0s2+  exp[(x  -  xr)lQ(x  -  xr)/(2m)]s2.  (15) 

xv 

The  first  term  under  the  radical  arises  from  the  error  in  the  estimate  (5  which  propagates  into 
yo  through  z^P,  while  the  second  term  comes  from  the  variance  of  the  error  at  the  point 

Xq- 


This  formulation  differs  somewhat  from  the  standard  prediction  confidence  interval, 

y0  =  yo  ±  t.osV  Zo(ZtZ)’1ZQS2  +  s2  ,  (16) 

/s 

due  to  the  fact  that  the  variance  of  3  is  (Z'V^Z )_1o2  instead  of  the  usual  (ZlZ)_1a2  and  the 
variance  of  £o  is  w0<T  =  exp[(x0-xr)tQ(x0-xr)/(2m)]a2  instead  of  a21.  The  prediction 
confidence  limit  of  equation  (15)  is  derived  in  Appendix  B. 

The  model  is  considered  valid  within  the  region  where  half  of  the  90%  confidence 
interval  is  less  than  some  value  (Kl)  determined  by  the  user.  The  magnitude  of  this  value 
can  be  a  reflection  of  the  risk  of  a  poor  product  that  the  user  is  willing  to  take.  This 
constraint  becomes: 

t.05 Vz}j(Zlv'lz)  ZqS2  +  exp[(x  -  xjtQix  -  xr)/(2m)]s2  <  Kj  (17) 

The  second  criteria  concerning  the  distance  form  the  data  is  a  rule  of  thumb.  In  general  it 
is  unwise  to  extrapolate  far  from  the  data  that  created  the  model.  Even  though  a  model  may 
appear  to  fit  the  data  well,  this  does  not  necessarily  mean  that  the  model  will  still  be  good 
far  from  the  data  from  which  it  was  created. 

This  constraint  reflects  the  fact  that  we  do  not  know  the  true  form  of  the  equation 
relating  the  control  variables  to  the  output.  If  the  form  of  the  equation  were  known  then 


1  We  chose  to  use  the  prediction  of  y  confidence  interval  instead  of  the  more  common  expected  value  of  y 

confidence  interval,  which  for  our  model  would  be  ^0  =  ^0  ~  *.05'  ZM^  ^  ^  Z0S  because  the 
prediction  confidence  interval  is  a  function  of  the  error  at  the  prediction  point,  and  we  are  using  the  error 
variance  to  "explain"  the  lack  of  fit  of  the  model  away  from  the  reference  point. 
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this  constraint  would  be  unnecessary.  To  quantify  the  feeling  that  we  should  not 
extrapolate  far  from  the  data,  the  generalized  multivariate  squared  distance  is  employed. 
The  multivariate  squared  distance  is: 


(x  -\iyt\x  -ji)  (i8) 

where  X  is  the  covariance  matrix  of  x.  This  is  the  multivariate  analog  of  the  univariate 
squared  distance: 


[(x-|i)/o]2  =  (x-[i)to'2(x-|a)  (19) 

These  equations  measure  distance  from  the  population  mean  in  terms  of  standard  deviations 
(Johnson  and  Wichem,  1988).  The  user  then  determines  the  maximum  (K2)  that  this 
squared  distance  can  be.  Again  the  choice  of  this  constant  can  be  a  reflection  of  the  risk 
that  the  user  is  willing  to  take.  If  x  were  multivariate  normally  distributed  then  multivariate 
squared  distance  would  have  a  chi-squared  distribution.  X  is  not  multivariate  normally 
distributed,  however  in  order  to  develop  a  measurement  of  how  far  away  is  too  far,  we 
chose  to  use  Z'Xq.io  to  define  the  limits  of  the  region  in  which  an  extrapolation  would  be 
considered  good.  Since  we  do  not  know  I  we  replace  it  in  equation  (18)  by  its  estimate  S . 

/s 

Likewise  is  substituted  for  |i.  We  arrive  at  the  second  constraint : 


(x  -S/S'^x  -4)  <  K2;  where  K2  =  2^0  l0.  (20) 

These  two  constraints  determine  the  region  in  which  the  model  of  the  process  is  considered 
valid. 

2.3  Mathematical  Statement  of  the  Problem 

We  may  now  state  the  problem  as: 


Maximize  y  =  zlP 


subject  to: 


1 05 V z^ZV'Z)  ^s2  +  exp[(x  -  xr)*Q(x  -  xr)/(2m)]s2  <  Kj  (21) 
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(x  -ite-Hx  -£)  <  2xo.io 

and  constraints  imposed  by  the  user,  such  as 
X2  >  3 

xi  +  X3  ^  7 


A  nonlinear  optimizer  solves  this  problem.  The  solution  is  the  next  point  recommended  by 
the  sequential  optimization  program.  The  nonlinear  optimization  algorithm  is  discussed  in 
Appendix  C. 
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3  APPLICATION  OF  SEQUENTIAL  OPTIMIZER  TO  PAPER 
HELICOPTER  DESIGN 


For  the  purpose  of  testing  the  value  of  using  physically  based  models  in  a 
sequential  optimization  algorithm,  the  sequential  optimizer  was  applied  to  the  design  of  a 
paper  helicopter. 

3,1  Description  of  System 

When  a  paper  helicopter  of  the  type  shown  in  Figure  2  is  dropped,  it  begins  to  spin  as 
it  falls  downward.  The  spinning  breaks  the  free-fall  of  the  helicopter  and  prolongs  the  time 
the  helicopter  remains  aloft.  The  three  design  parameters  which  affect  flight  time  are  the 
length  of  the  blade  (L),  the  width  of  the  wing  (b),  and  the  total  weight  (wt)  of  the 
helicopter.  The  total  weight  is  equal  to  the  weight  of  the  paper  (wp)  plus  the  weight  added 
(wa)  at  the  end  of  the  tail.  The  angle  at  which  the  blades  are  initially  bent  was  determined 
to  have  negligible  affect  during  early  screening  experiments.  We  chose  to  design  the 
helicopters  for  long  flight  time  or,  equivalently,  low  velocity  (V).  (Actually,  negative 
velocity  was  used  as  the  objective  function  since  the  simple  optimizer  was  set  up  for 
maximization.) 


Figure  2:  Helicopter  design 
3.2  Physical  Model 

Dimensional  analysis  is  a  technique  commonly  used  in  the  physical  sciences  to 
simplify  the  analysis  of  complex  multivariable  problems  by  grouping  variables  into 
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dimensionless  quantities.  The  underlying  basis  for  dimensional  analysis  is  that  physical 
results  must  be  independent  of  the  system  of  units  chosen  for  their  measurement. 

The  Buckingham  Pi  theorem  (Bridgman,  1931)  provides  a  step  by  step  approach  to 
the  formation  of  sets  of  dimensionless  quantities  which  are  appropriate  to  a  specific 
problem.  The  set  of  dimensionless  variables  is  smaller  in  number  that  the  original  set  of 
variables.  In  the  case  of  the  helicopter,  there  are  four  variables  which  influence  the  velocity 
of  the  helicopter;  L,  b,  w,  and  p,  where  p  is  the  density  of  the  air.  We  may  formulate  an 
implicit  relationship  between  V  and  the  four  variables  as  follows: 


g(V,  L,  b,  w  ,  p)  =  0  (22) 

An  application  of  the  Pi  Theorem  results  in  the  transformation  of  these  five 
variables  into  two  dimensionless  groupings  Vv^b/VW  ancj  L/b,  from  which  an  implicit 
relationship  may  be  formed  which  characterizes  the  helicopter  problem: 

g(VVp  b/Vw,  L/b)  =  0  (23) 

This  relationship  may  be  restated  as: 


VVp  b/Vw  =  f(L/b) 


(24) 


We  chose  to  approximate  the  function  of  (L/b)  with 

f(L/b)  =  p0  +  PjL/b  +  p2(L/b)2  (25) 


so  that 


Vfp  b/Vw  =  Po  +  Pi  L/b  +  P2(L/b)2.  (26) 

Since  p  is  constant  it  will  be  absorbed  into  the  coefficients.  The  model  then  becomes 

V  =  Vw/b[  po  +  PiL/b  +  P2(L/b)2].  (27) 

The  application  of  dimensional  analysis  to  the  helicopter  problem  has  provided 
benefit  in  two  ways.  First,  the  number  of  experimental  control  factors  has  been  reduced 
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from  the  original  list  of  L,b,  and  w  to  VtfptyVw  and  L/b.  Second,  the  grouping  of  the 
variables  into  the  dimensionless  terms  and  the  grouping  of  the  teims  themselves  stems  from 
physical  considerations,  and  may  therefore  be  expected  to  model  the  process  more 
accurately  than  a  simple  polynomial  model. 

A  by-product  of  the  fact  that  the  number  of  experimental  parameters  has  been 
reduced  from  3  to  2  is  that  specification  of  the  two  dimensionless  parameters  Wp’b/Vw  and 
IJb  does  not  uniquely  specify  the  three  helicopter  design  parameters  L,  b,  and  w.  Indeed, 
if  the  model  resulting  from  dimensional  analysis  captured  all  physical  effects  (which  it  does 
not),  the  two  dimensionless  parameters  would  be  sufficient  to  predict  the  performance  of 
the  design.  Since  the  the  dimensional  analysis  is  an  approximate  physical  model,  an 
additional  term,  L,  is  added  to  the  model  so  that  L,  b,  and  w  are  independently  determined. 
The  final  physically  based  model  is: 


V  =  Vw/b[  Po  +  PiUb  +  p2(L/b)2]  +  p3L. 


(28) 


3.2  Application  of  Basic  Sequential  Optimizer 

The  sequential  optimizer  was  used  to  improve  the  design  of  the  helicopters.  The 
optimization  process  was  performed  with  the  three  different  models  shown  below: 


Model  1:  V  =  Vw/b[  p0  +  PiUb  +  p2(L/b)2]  +  p3L 


(29) 


Model  2:  V=Po  +  PiL  +  P2b  +  P3w 


(30) 


Model  3:  V  =  Po  +  PiL  +  p2b  +  p3w  +  P4L2  +  p5b2  +  Pgw2+  p2Lb+  PgLw+P9bw  (31) 

The  first  model  is  the  physical  model  based  on  the  dimensional  analysis.  The  second 
model  is  a  model  which  is  linear  in  the  basic  variables  L,  b,  and  w.  Model  1  and  Model  2 
both  have  four  terms.  The  third  model  is  a  full  quadratic  in  the  basic  variables. 


In  order  to  start  the  optimization  process,  data  sufficient  to  estimate  the  coefficients  and  the 
error  variance  of  the  models  were  required.  The  starting  data  were  obtained  through 
parallel  DOE.  The  sequential  optimizer  was  nin  using  two  separate  sets  of  starting  data. 
Set  1,  and  Set  2,  shown  in  Table  1  and  Table  2,  respectively.  This  wa:  done  so  that  the 
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effect  of  the  models  on  the  optimization  speed  could  be  compared  under  different  starting 
conditions.  Model  1  and  Model  2  were  compared  based  on  the  first  set  of  starting  data. 
The  first  set  of  starting  data  came  from  a  two  level  full  factorial  DOE. 


RUN# 

L 

Set  1 

b 

wt 

V  ave 

(cm) 

(cm) 

(g) 

(cm/s) 

1 

8.0 

2.5 

10.0 

338.0 

2 

12.0 

2.5 

10.0 

187.8 

3 

8.0 

3.5 

10.0 

265.5 

4 

12.0 

3.5 

10.0 

181.9 

5 

8.0 

2.5 

12.0 

390.3 

6 

12.0 

2.5 

12.0 

201.8 

7 

8.0 

3.5 

12.0 

313.8 

8 

12.0 

3.5 

12.0 

184.0 

Table  1 :  First  set  of  starting  data 


All  three  models  were  then  compared  based  on  the  second  set  of  starting  data.  Since 
there  were  initially  only  six  data  points  the  full  third  model  could  not  be  estimated.  Due  to 
this  lack  of  data  points,  only  the  first  four  terms  of  the  model  (coefficients  0  -  3)  were  used 
at  the  first  step,  and  an  additional  term  from  the  full  quadratic  model  was  added  after  each 
additional  experiment. 


Set  2 


RUN# 

L 

b 

wt 

V  ave 

(cm) 

(cm) 

(g) 

(cm/s) 

1 

9.0 

2.0 

10.0 

344.9 

2 

7.0 

3.0 

10.0 

329.8 

3 

9.0 

3.0 

12.0 

290.0 

4 

8.0 

2.5 

11.0 

356.5 

5 

8.0 

3.0 

11.0 

303.4 

6 

9.0 

3.0 

10.0 

260.3 

Table  2:  Second  set  of  starting  data 


For  this  problem  there  were  two  user  constraints  on  the  control  variables.  Since  the 
paper  from  which  the  helicopters  were  constructed  was  25  cm  long,  the  length  of  the 
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helicopter,  excluding  the  tail,  was  constrained  to  be  less  than  25  cm.  Also,  the  lower 
bound  on  helicopter  weight  was  set  at  8  grams.  This  weight  bound  was  chosen  so  that  the 
optimizer  would  never  recommend  a  physically  impossible  design  -  one  where  the  weight 
of  the  paper,  wp,  in  the  helicopter  which  is  determined  by  L  and  b  would  be  greater  than 
the  recommended  weight,  wt  The  constrained  optimization  problem  was 

Maximize  (-V) 

subject  to:  (32) 

L  +  b  <  25 

w  >  8 

The  user  specified  parameters  were  set  as  follows.  The  A;,  the  scaling  constants  of 
the  error  weighting  function  of  equation  (4), were  chosen  to  be  the  range  of  their  respective 
control  variables  x;  in  the  initial  data  set.  This  means  that  for  the  the  experiments  which 
began  with  the  first  initial  data  set,  Ai  =  4,  A2  «  1,  and  A3  =  2.  For  the  second  initial  data 
set  Ai  =  2,  A2  =  1,  A3  =  2.  Kl,  the  variable  which  determines  the  maximum  error  in  the 
prediction  of  y  as  specified  in  equation  (17),  is  also  set  by  the  user.  For  the  design 
processes  based  on  the  first  set  of  initial  data,  KI  was  set  at  67.1  for  the  entire  run  of  the 
optimization.  For  the  design  processes  based  on  the  second  sei  of  initial  data,  Kl  was 
originally  set  at  23.5  cm/s.  However,  after  the  first  2  sequential  experiments,  Kl  had  to  be 
changed  to  67.1  cm/s  due  to  the  high  predictive  uncertainty  in  Models  2  and  3.  K2,  the 
bound  for  the  distance  constraint,  was  set  at  2x2oi  =  12.5  for  all  of  the  optimizer 
applications. 

Data  on  the  descent  velocity  were  obtained  by  releasing  the  helicopters  in  a  stairwell 
from  a  height  of  three  stories  and  the  time  of  flight  was  used  to  calculate  velocity.  Because 
of  variability  in  the  flight  time,  replicates  of  each  experiment  were  performed  and  averaged 
to  become  one  data  point  in  the  optimization.  The  "velocity"  entered  into  the  sequential 
optimizer  was  calculated  as  100/time,  however  the  actual  velocity  (in  terms  of  cm/s)  is 
recorded  in  the  tables  of  this  thesis.  An  occasional  clearly  outlying  point  was  not 
considered  in  the  determining  the  average  velocity.  After  the  starting  data  was  collected, 
the  sequential  optimizer  was  then  used  to  determine  the  subsequent  helicopter  designs. 
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4  RESULTS  AND  DISCUSSION 


4.1  Experimental  Results 

The  optimizer  was  run  with  five  model  and  initial  data  set  combinations,  specifically:. 
Model  1  and  Set  1,  Model  2  and  Set  1,  Model  1  and  Set  2,  Model  2  and  Set  2,  and  Model  3 
and  Set  2.  The  presentation  of  the  data  is  organized  according  to  the  initial  data  set 

The  results  of  the  comparison  of  Model  1  and  Model  2  based  on  the  Set  1  of  starting 
data  are  presented  in  Figure  3.  The  first  eight  data  points  are  from  the  initial  data  set 
obtained  through  parallel  design,  and  the  remaining  data  points  are  the  result  of  the 
sequential  optimizations. 


Sequential  Optimization  Based  on  Set  1  of  Initial  Data 


Experiment  # 


Figure  3:  Optimizations  based  on  the  first  set  of  initial  data 

As  may  be  seen  by  referring  to  Figure  3,  the  performance  of  the  optimizer  with  the 
dimensional  analysis  model  (Model  1 )  was  found  to  be  superior  to  the  performance  of  the 
optimizer  with  the  linear  model.  The  optimizer  based  on  the  dimensional  analysis  model 
reached  a  point  near  the  optimum  by  the  fourth  sequential  design.  At  this  point  the  next 
recommended  design  differed  only  slightly  from  the  last  and  we  chose  to  treat  this  design 
as  the  optimum.  The  optimization  based  on  the  linear  model  had  not  yet  reached  the 
optimum  design  by  the  tenth  sequential  design,  although  it  had  also  improved  the  helicopter 
design. 
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The  step  by  step  evolution  of  the  helicopter  designs  in  L-b  design  space  illustrates  the 
way  in  which  the  sequential  optimizer  functions.  The  location  of  the  designs  for  the 
optimizations  based  on  Set  1  are  shown  in  Figures  4  and  5.  The  shaded  region  represents 
the  infeasible  designs,  where  L  +  b  >  25. 


Figure  4:  Design  locations  for  the  optimization  based  on  Model  1. 


Figure  5:  Design  locations  for  the  optimization  based  on  Model  2. 
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The  optimum  design  appears  to  be  on  the  bounds  of  the  design  region  -  where  L  +  b  =25, 
and  w  =  8. 

The  actual  values  of  the  data  plotted  in  Figures  4  and  5  are  listed  in  Tables  3  and  4. 
Additionally,  the  predicted  value  of  velocity  ,  as  well  as  half  of  the  width  of  the  90% 
confidence  interval  for  the  prediction  are  given.  Note  that  the  prediction  error  constraint  of 
equation  (17)  requires  that  half  of  the  width  of  the  90%  confidence  interval  be  less  than  Kl, 
which  is  67.1  cm/s  for  the  design  processes  based  on  the  first  set  of  initial  data.  If  the 
optimum  for  the  step  was  on  the  bound  of  the  region  where  the  model  is  considered  valid, 
then  the  constraint  which  caused  the  bound  is  listed  in  the  constr  column.  The  letter  p 
indicates  that  the  prediction  error  constraint,  which  is  defined  by  equation  (17)  and  the 
choice  of  Kl,  was  the  binding  constraint,  while  d  indicates  that  the  distance  constraint, 
defined  by  equation  (20)  and  the  choice  of  K2,  was  the  binding  constraint. 


Physical  Dimensional  Anaysis  Model 


L 

b 

wt 

con¬ 

V  ave 

V  pred 

± 

(cm) 

(cm) 

(g) 

straint 

(cm/s) 

(cm/s) 

(cm/s) 

17.4 

3.0 

9.9 

d,p 

157.7 

13.2 

67.1 

18.0 

3.3 

8.0 

p 

120.0 

90.3 

67.1 

20.6 

4.4 

8.0 

r 

d 

116.8 

65.4 

66.8 

20.7 

4.3 

8.0 

★ 

110.6 

109.2 

29.5 

Table  3:  Results  of  the  sequential  optimization  based  on  Model  1  and  the  first  set  of  initial 
data 
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Linear  Model 


RUN  # 

L 

b 

wt 

con¬ 

V  ave 

V  pred 

± 

(cm) 

(cm) 

(g) 

straint 

(cm/s) 

(cm/s) 

(cm/s) 

1 

13.7 

3.4 

10.7 

P 

171.6 

113.9 

67.1 

2 

13.8 

3.3 

10.5 

P 

171.6 

137.3 

67.1 

3 

15.0 

3.4 

10.4 

P 

166.2 

110.8 

67.1 

4 

15.5 

3.4 

10.4 

P 

155.6 

118.1 

67.1 

5 

17.0 

3.6 

10.2 

P 

138.9 

93.1 

67.1 

6 

19.5 

3.8 

10.0 

P 

120.7 

61.8 

67.1 

7 

21.4 

3.6 

8.7 

P 

109.2 

65.4 

67.1 

8 

20.7 

4.3 

11.6 

P 

123.2 

80.3 

67.1 

9 

21.8 

3.2 

8.0 

P 

138.8 

83.4 

67.1 

10 

20.2 

4.8 

11.4 

P 

107.8 

91.1 

67.1 

Table  4:  Results  of  the  sequential  optimization  based  on  the  Model  2  and  the  first  set  of 
initial  data 


Run  eight  was  repeated  for  the  linear  model  and  the  second  data  point  was  used  in  the 
sequential  optimization  process. 


The  next  group  of  results  comes  from  the  applications  of  the  sequential  optimizer 
based  on  the  second  set  of  starting  data.  All  three  models  were  compared  based  on  this  set 
of  starting  data.  Figure  6  shows  the  results  of  all  three  design  processes.  The  first  six  data 
points  are  from  the  initial  data  set. 

Sequential  Optimization  Based  on  Set  2  of  Initial  Data 


Q  Initial  Data 
— Model  1 
•  Model  2 
Model  3 


Figure  6:  Optimizations  based  on  the  second  set  of  initial  data 
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The  comparison  of  Model  1  and  Model  2  based  on  the  second  set  of  initial  data  gives 
results  similar  to  those  based  on  the  first  initial  data  set.  The  physically  based  model 
improves  the  design  the  most,  although  both  models  appears  to  improve  the  helicopter 
design. The  third  model,  which  adds  a  new  term  from  the  full  quadratic  at  each  step,  gave 
the  weakest  performance,  although  it  too  appeared  to  improve  the  design.  None  of  the 
models  reached  a  final  design  by  the  end  of  the  six  sequential  designs. 


The  numerical  results  of  the  three  design  processes  are  listed  in  tables  4-6. 

Physical  Dimensional  Analysis  Model 


RUN  # 

L 

b 

wt 

(g> 

con¬ 

V  ave 

V  pred 

± 

(cm) 

(cm) 

straint 

(cm/s) 

(cm/s) 

(cm/s) 

1 

8.6 

3.0 

8.5 

P 

225.7 

293.3 

23.5 

2 

8.9 

3.1 

8.0 

P 

219.3 

204.9 

23.5 

3 

10.2 

3.6 

8.0 

d 

204.9 

149.4 

30.3 

4 

11.4 

4.0 

8.0 

d 

162.8 

157.9 

48.5 

5 

12.8 

4.5 

8.0 

d 

151.7 

130.5 

35.0 

6 

14.6 

5.1 

8.0 

d 

129.8 

121.1 

25.6 

Table  5:  Results  of  the  sequential  optimization  based  on  the  Model  1  and  the  second  set  of 
initial  data 


Linear  Model 


RUN# 

L 

b 

wt 

(g) 

con¬ 

V  ave 

V  pred 

+ 

(cm) 

(cm) 

straint 

(cm/s) 

(cm/s) 

(cm/s) 

1 

9.3 

3.5 

10.7 

P 

225.3 

214.7 

23.5 

2 

10.2 

4.2 

10.9 

d 

212.7 

137.3 

22.5 

3 

11.1 

4.8 

10.9 

P 

212.3 

141.5 

67.1 

4 

11.9 

5.4 

10.9 

P 

118.3 

158.2 

67.1 

5 

14.1 

7.1 

11.1 

d 

169.7 

107.8 

66.3 

6 

16.4 

8.2 

10.7 

d,p 

167.8 

119.7 

67.0 

Table  6:  Results  of  the  sequential  optimization  based  on  the  Model  2  and  the  second  set  of 
initial  data 
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Quadratic  Model 


RUN# 

L 

b 

wt 

con¬ 

V  ave 

V  pred 

± 

(cm) 

(cm) 

(g) 

straint 

(cm/s) 

(cm/s) 

(cm/s) 

1 

9.3 

3.5 

10.7 

P 

225.3 

214.7 

23.5 

2 

9.3 

3.7 

10.7 

P 

249.7 

204.9 

23.5 

3 

9.6 

3.7 

9.6 

P 

235.2 

206.3 

67.1 

4 

9.8 

3.6 

9.8 

P 

224.6 

191.9 

67.1 

5 

10.0 

3.4 

9.8 

P 

218.1 

203.0 

67.1 

6 

10.3 

3.3 

10.0 

P 

221.6 

200.5 

67.1 

Table  7:  Results  of  the  sequential  optimization  based  on  the  Model  3  and  the  second  set  of 
initial  data 


The  distance  constraint  frequently  limited  the  physical  model  while  the  polynomial  models 
were  limited  by  the  constraint  that  the  prediction  confidence  interval  be  less  than  23.5  for 
the  first  two  sequential  designs  and  less  than  67.1  for  the  remaining  sequential  designs. 
Run  number  5  was  repeated  for  both  the  linear  and  quadratic  models.  The  sequential 
optimization  process  proceeded  on  the  basis  of  the  repeated  runs. 

4.2  Discussion  of  Model  Comparison 

The  physically  based  dimensional  analysis  model  performed  the  best  for  both  sets  of 
initial  data.  The  improvement  in  the  helicopter  design  after  a  given  number  of  steps  was 
greater  for  the  optimization  based  on  the  dimensional  analysis  model  than  it  was  for  the 
optimization  based  on  either  of  the  polynomial  models.  This  can  be  attributed  to  two 
interconnected  effects  which  reflect  the  fact  that  the  physical  model  is  a  better  representation 
of  the  true  system.  The  first  is  that  the  physically  based  model  provided  a  better  fit  to  the 
data  so  that  the  region  in  which  the  model  is  considered  valid  was  larger.  The  second  is 
that  the  model  provides  for  better  direction  within  the  region.  This  region,  where  both  the 
distance  constraint  and  the  prediction  error  constraint  are  satisfied,  is  larger,  because  the 
prediction  error  constraint  is  less  limiting.  (The  distance  constraint,  which  is  a  function  of 
the  distribution  of  the  control  variables,  is  independent  of  the  model,  while  the  prediction 
error  constraint  depends  on  how  well  the  model  fits  the  data.)  Extrapolations  based  on  the 
model  are  believed  to  be  sufficiently  accurate  within  the  region  where  the  model  is  valid. 
The  sequential  optimizer  does  recommend  designs  beyond  this  point.  Because  of  the  larger 
region,  the  design  could  change  more  between  each  step,  speeding  the  sequential 
optimization  process. 
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The  physical  model  frequently  placed  the  step  optimum  on  the  bound  created  by  the 
distance  constraint,  not  the  prediction  error,  indicating  that  poor  fit  of  the  model  did  not 
limit  the  optimization.  The  step  optimum  was  on  the  distance  bound  for  6  of  10  steps  for 
the  physical  model.  (The  step  optimum  was  at  the  intersection  of  the  distance  and 
prediction  error  constraints  at  one  of  these  points.)  This  was  in  part  due  to  the  good  fit  of 
the  physical  model  to  the  data.  Another  factor  is  that  if  the  optimum  was  at  the  distance 
bound  at  the  previous  step  it  may  have  the  propensity  to  be  there  again  since  the  distance 
constraint  may  not  change  significantly  between  steps,  with  the  addition  of  one  data  point. 
The  distance  is  measured  in  standard  deviations  of  the  data,  and  the  standard  deviation  may 
not  change  much  with  the  addition  of  one  point.  The  polynomial  models,  on  the  other 
hand,  rarely  reached  the  distance  constraint,  being  limited  by  the  prediction  error  of  the 
model.  The  step  optimum  for  the  polynomial  models  was  on  the  distance  bound  only  3  of 
22  times  and  one  of  these  points  was  at  the  intersection  of  the  prediction  error  ond  the 
distance  constraints. 

Although  the  simple  linear  model  performs  well  for  the  sequential  optimization  of  the 
paper  helicopter,  linear  models  in  general  will  not  be  able  to  optimize  processes  or  designs 
where  the  relationship  betw-een  the  response  and  the  basic  variables  contains  curvature.  A 
quadratic  model  would  resolve  this  problem. 

The  quadratic  model,  which  contains  the  linear  model  as  a  subset  of  its  terms, 
performs  poorly  due  to  its  high  predictive  error,  which  restricts  the  region  in  which  the 
model  is  valid.  This  is  due  to  the  large  number  of  estimated  coefficients.  Each  estimated 
coefficient  has  a  variance  as  well  as  a  covariance  with  the  each  of  the  other  estimated 
coefficients.  The  variance  and  covariance  of  each  estimated  coefficient  propagates  into  the 
variance  of  y. 

4.3  Suggested  Improvements  to  the  Sequential  Optimization  Algorithm 

The  optimization  algorithm  performed  well,  improving  the  design  of  the  paper 
helicopter  for  both  starting  data  sets  and  all  three  models.  However,  the  algorithm  can  be 
improved  in  a  number  of  areas. 

One  aspect  of  the  current  algorithm  which  is  unsatisfactory  is  that  if  the  model  is  a 
poor  fit  to  the  data  and  the  prediction  error  is  high,  then  there  may  be  not  be  a  region  where 
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the  model  is  valid.  The  optimizer  cannot  function  in  this  situation  unless  the  bound  on  the 
maximum  prediction  error  (Kl)  is  increased,  meaning  that  less  accurate  predictions  are 
acceptable.  This  happened  during  the  optimization  processes  based  on  Model  2  and  3  and 
the  starting  data  from  Set  2.  Kl  was  originally  set  at  23.5,  but  then  had  to  be  changed  to  a 
larger  value  after  two  sequential  designs.  Though  this  was  not  critical  in  the  helicopter 
application  of  this  paper,  in  many  on-line  applications  where  avoiding  the  production  of 
scrap  is  a  primary  concern,  the  accuracy  of  the  model  predictions  is  important.  One 
solution,  exemplified  by  the  physically  based  model,  is  the  usage  of  better  models. 

Another  solution  for  the  case  where  a  better  model  is  not  available,  is  the  shrinkage  of 
the  local  region,  so  that  the  region  where  a  good  fit  to  the  data  is  required  is  smaller,  since 
the  problem  results  from  the  model  being  insufficient  to  explain  the  data  over  the  entire 
weighted  range.  Linear  and  quadratic  models,  in  particular  will  be  accurate  over  a 
sufficiently  small  region.  This  could  be  accomplished  by  decreasing  one  or  all  of  the  Aj's . 
However,  this  may  result  in  models  created  from  few’  effective  points.  This  problem  could 
be  resolved  by  adding  an  additional  constraint  that  the  model  be  based  on  a  minimum 
number  of  effective  data  points.  When  this  constraint  is  violated,  the  number  of  effective 
points  could  be  increased  by  doing  experiments  which  are  close  to  the  reference  point. 

Additional  suggestions  for  improving  the  algorithm  include  using  the  best  data  point 
instead  of  the  last  data  point  as  the  reference  point  around  which  the  model  is  created.  A 
further  improvement  would  be  the  development  of  a  starting  algorithm  which  would 
eliminate  the  need  for  a  parallel  design  of  starting  data.  The  basic  idea  of  a  starting 
algorithm  would  be  to  begin  optimization  process  while  the  initial  set  of  data  is  being 
collected. 
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5  CONCLUSIONS 


In  this  work  we  have  sought  to  explore  the  benefits  of  incorporating  physically 
based  models  into  a  sequential  optimization  algorithm.  We  hypothesized  that  physically 
based  models  which  incorporate  an  engineer's  prior  knowledge  of  a  system  would  lead  to 
more  rapid  optimization  of  products  and  processes. 

In  order  to  test  this  idea,  we  developed  a  simple  sequential  optimizer  which  could 
be  modified  to  : a  ccommodate  different  models.  The  optimization  algorithm  is  based  on  the 
sequential  use  of  local  models  created  through  weighted  linear  regression.  A  new  operating 
point  or  design  is  recommended  at  the  optimum  of  the  model  within  the  region  where  the 
model  is  considered  valid.  Extrapolations  based  on  the  model  are  believed  to  be  accurate 
w  ithin  this  region.  A  new  model  is  created  after  each  new  data  point  is  collected. 

As  a  test  case,  we  decided  to  design  a  paper  helicopter  for  maximum  flight  time. 
The  optimizer  was  run  using  a  physically  based  model  and  simple  polynomial  models  in 
order  to  compare  them.  The  physically  based  model  was  derived  by  dimensional  analysis, 
a  technique  that  groups  variables  based  on  their  dimensions.  The  use  of  the  physically 
based  model  derived  from  dimensional  analysis  significantly  reduced  the  number  of  steps 
required  to  optimize  the  helicopter  design. 

Our  work  indicates  that  the  use  of  physically  based  models  in  the  context  of 
sequential  optimization  leads  to  the  rapid  optimization  of  processes  and  product  designs. 
We  feel  that  this  approach  is  especially  useful  for  the  optimization  of  manufacturing 
processes.  Future  work  might  include  the  development  and  verification  of  physically  based 
models  for  important  classes  of  processes.  One  can  envision  equipment  controllers  that 
have  embedded  in  them  optimization  capabilities,  each  with  a  process  specific  model.  The 
capability  to  rapidly  optimize  processes  will  become  more  valuable  as  the  trend  toward 
flexible  manufacturing  with  its  shorter  runs  and  change  over  times  continues. 
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APPENDIX  A 


DERIVATION  OF  WEIGHTED  REGRESSION 

A  brief  derivation  of  the  formula  for  weighted  regression  is  given  in  this  appendix. 
For  a  more  thorough  discussion  see  Draper  and  Smith. 

Weighted  regression  is  used  when  the  error  at  each  observation  does  not  (or  is 
assumed  to  not)  conform  to  the  standard  linear  regression  assumption  that  the  error 
variance  of  each  observation  is  equal  and  independent.  Weighted  regression  basically 
consists  of  transforming  the  variables  involved  in  the  regression  so  that  these  assumptions 
are  satisfied.  This  is  done  by  premultiplying  the  standard  linear  regression  model  by  a 
matrix  so  that  the  resulting  error  conforms  to  the  standard  regression  assumptions. 

In  the  context  of  least  squared  regression  it  amounts  to  penalizing  the  squared  error 
differently  for  different  points.  The  model  of  the  situation  where  weighted  regression  is 
used  is: 

y  =  zp  +  C  (A.i) 

where 

C  ~  N(0,Va2)  (A.2) 

Since  V  is  square  and  symmetric  P  can  calculated  such  that 

ptp  =  Pp  =  P2  =  v.  (A. 3) 

Premultiplying  the  linear  model  by  P'1  yields: 

Ply  =  P!zp  +P-1;.  (A.4) 

Let  the  following  variables  be  defined  by 

y*  =  P^y,  Z*  =  P^Z,  and  C*  =  P1;  (A.5) 

Thus  equation  (A.4)  can  be  rewritten  as 
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y*  =  z*p  +  C  . 


(A. 6) 


where  the  distribution  of  C  is 


C*  ~  N(0,lc2) 

(This  is  because  Cov(Cx)  =  CCov(x)Ct  when  C  is  a  constant  matrix  and 
variable  vector.)  The  standard  regression  formula  can  now  be  applied, 
regression  equation  is 


(3=  (Z‘Z)-1Z*y 

By  substituting  the  starred  quantities  into  this  equation  we  arrive  at 

(3=  (z-'zv'zV 

Returning  to  the  original  variables,  the  equation  is 

P=  [(P^Z/fP^Z)]  1(P*1Z)t(P‘1y). 
Since  (AB)1  =  BlAl  the  expression  can  be  rewritten  as 

P=  [ZKP^kP^Z)]'^^1)^-^) 

After  applying  the  fact  that  (A1)'1  =  (A'1)1  the  equation  becomes 

P=  [ZW^Z^ZW^y) 


Since  Pl  =  P 


p=  [Z'P^P^Z^^Z'P-HP’V) 


(A.7) 

x  is  a  random 
The  standard 


(A. 8) 


(A. 9) 


(A.  10) 


(A.  11) 


(A.  12) 


(A.13) 
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which  reduces  to 


p=  (Z'V^Zr’Z'V'V  (A. 14) 


The  covariance  of  (3  is  required  for  the  calculation  of  confidence  intervals.  The  covariance 
is  derived  as  follows: 


Cov(p)  =  Cov[(ZtV-1Z)'1ZtV4y]  (A.  15) 

The  first  step  is  applying  the  relationship  Cov(Cx)  =  CCov(x)Cl  to  equation  (A.  15).  This 
results  in 


A 

Cov(p)  =  (Z,V‘1Z)‘1ZtV'1Cov(y)[(ZrV"1Z)'1ZtV'1]t  (A.16) 

The  expression  of  equation  (A.16)  can  be  simplified.  The  next  3  steps  are  basic  matrix 
manipulation  applying  (AB)1  =  B'A£  and  (A1)’1  =  (A4)1  to  equation  (A.16). 


Cov(P)  =  (ZrV'1Z)"1ZtV'1Cov(y)(V‘1)tZ[(ZtV4Z)t]  (A.17a) 

Cov(p)  =  (ZtV1Z)"1ZlV'1Cov(y)(Vt)‘1Z[Zt(V'1)lZ]  1  (A.17b) 

Cov(p)  =  (ZtV'1Z)’1ZtV'1Cov(y)(Vt)'1Z[Zt(Vt)’1Z]  1  (A.  17c) 

finally,  since  V1  =  V,  the  equation  becomes: 

Cov(p)  =  (ZlV'1Z)'1ZtV'1Cov(y)V'1Z(ZtV'1Z)'1  (A.  18) 

Substituting  y  =  ZP  +  ^  into  equation  (A.  18)  yields 

A 

Cov(p)  =  (ZtV'1Z)'1ZtV'1Cov(ZP  +  OV^ZfZ'V^Z)'1  (A.19) 

Since  ZP  is  constant 


33 


• 

Cov(p)  =  (ZtV-1Z)‘1ZtV-1Cov(OV-1Z(ZtV-1Z)'1. 

(A. 20) 

Substituting  Cov(^)  =  Vo2  into  equation  (A. 20)  yields 

Cov(p)  =  (ZtV-1Z)‘1ZtV1Vo2V-1Z(ZtV1Z)'1 

(A. 21) 

which  can  be  rearranged  as 

Cov(p)  =  (ZtV'1Z)‘1ZtV'1Z(ZtV'1Z)"1a2 

(A. 22) 

and  simplified  to 

Cov(P)  =  ylo2. 

(A. 23) 

XV 

Alternately,  this  derivation  can  be  arrived  at  by  noting  that  the  covariance  of  P 
unweighted  regression  is: 

in  an 

• 

Cov(p)  =  (Z'zy1©2. 

(A. 24) 

Substituting  in  the  the  starred  variables  yields: 

Cov(p)  =  (z^zyV 

(A.25) 

Returning  to  the  original  variables  the  equation  becomes 

Cov(P)  =  [(P^Z/P-'Z  ]lo2 

(A.26) 

which  after  the  following  matrix  manipulations 

Cov(P)  =  [(ZP-VP^ZfV 

(A. 27 a) 

Cov(p)  =  [Z'IP-Vp-’ZJ'V 

(A.  27b) 

• 

Cov(p)  =  (Z'P^P'Z)'1©2 

(A  .27c) 
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becomes 


Cov((3)  =  (Z'V-1Z  )'lc2.  (A. 28) 

The  covariance  for  P  from  an  unweighted  regression  is  derived  in  Johnson  and  Wichem 

and  the  covariance  of  (3  for  a  weighted  regression  is  given  in  Draper  and  Smith. 

The  variance,  a2,  is  estimated  by 

t  /s 

s2  =  (y-Xp)V1(y-Xp)/(n-p)  (A.29) 

where  n  is  the  number  of  data  points  and  p  is  the  number  of  estimated  coefficients  in  the 
model. 
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APPENDIX  B 


DERIVATION  OF  PREDICTION  CONFIDENCE  INTERVAL  FOR  WEIGHTED 
REGRESSION 

The  prediction  confidence  interval  for  the  linear  model  of  this  thesis,  is  derived  in  this 
appendix.  The  prediction  confidence  interval  for  the  case  where  P  is  estimated  through 
standard  unweighted  regression  and  error  is  assumed  uniform  at  all  points  is  derived  in 
Johnson  and  Wichem.  According  to  the  model  the  response  at  the  point  zo  is 


>'o  =  zbP  +  £o 


(B.l) 


which  is  estimated  by 


y0  =  zoP- 


(B.2) 


Thus  the  error  in  the  prediction  is 


yo  -  yo  =  (zbP  +  £o  -  z^P)  (B.3) 

We  observe  that  since  P  is  normally  distributed,  the  linear  combination  zbP  is  also 
normally  distributed.  Further,  eq  has  a  normal  distribution,  so  the  quantity  y0  -  yo  i  s 
normally  distributed.  Now  we  must  find  the  mean  and  the  variance  of  y0-  yo.  The 

expected  value  of  y0  -  yo  is 


E(y0  -  yo)  =  E(zbP  +  eo  -  ^P)  =  ^P  +  0  -  zbP  =  0  (B.4) 


Since  ZqP  is  constant  and  zbp  and  Eo  are  independent 

var(y0  -  yo)  =  var(z‘0p  +  Eo  -  z^P)  =  var(eo)  +var(zbp)  (B.5) 


At  this  point  the  derivation  deviates  from  the  standard  derivation  for  an  unweighted 
regression.  The  variance  of  £o  is 
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var(e0)  =  w0a2  =  exp[(x0-xr)tQ(x0-xr)/(2m)]a2 


(B.6) 


and  the  variance  of  zoP  is 

/V  /V 

varfz^P)  =  zJ)Cov(p)zo  =  z^fZ'V^Z)'1^2  (B.7) 

Substituting  equations  (B.6)  and  (B.7)  into  (B.5),  the  variance  of  y0  -  9o  can  be  written 
as 


var(yo-yo)  =  z^Z'V^ZJ^zoO2  +  w0o2  (B.8) 

The  first  term  represents  the  variance  due  to  the  error  at  point  xo,  and  the  second  term 
represents  the  variance  due  to  the  error  in  estimating  P  We  can  now  write  the  distribution 

of  y0  *  9o  as 


yo  -  yo  ~  N(0,Zo(ZtV‘1Z)'1zoo2  +  w0o2)  (B.9) 

From  this  distribution  of  y,  it  can  be  shown  that 

(yo-yo)A/,z))(ZtV-1Z)'1zos2  +  w0s2  ~  t  (B .  1 0) 

The  confidence  interval  for  the  prediction  value  of  yo  is  then  derived  in  the  usual  manner 
yielding: 

yo  =  yo  ±  tot/2, V  z[)(ZtV'1Z)'1zos2  +  w0s2  (B.ll) 


or 


y0  =  y0  ±  ta/2"V  Zq(Z1V_1Z)  1z0s2+exp[(x0  -  xr)lQ(x0  -  xr)/(2m)]s2  (B.12) 

The  prediction  confidence  interval  for  the  case  where  P  is  estimated  through  unweighted 
regression  is  derived  in  Johnson  and  Wichem. 
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appendix  c 


EASIC  OPTIMIZATION  ALGORITHM 

The  optimization  problem  considered  in  this  thesis  is  of  the  form 

max  f(x) 

subject  to 

g(x)  >  c 
h(x)  >  b 

where 


(C.l) 


h(x)  is  the  set  of  linear  constraints. 
g(x)  is  the  set  of  nonlinear  constraints. 

In  our  algorithm  there  are  only  two  nonlinear  constraints.  They  which  arise  from  the 
constraints  which  determine  the  feasible  region,  since  the  constraints  that  the  user  imposes 
are  required  to  be  linear. 

The  algorithm  used  to  solve  this  problem  is  based  on  the  gradient  projection  method.  The 
algorithm  is  moves  from  the  current  feasible  point  xk  to  an  improved  feasible  point  xk+i. 
The  gradient  projection  method  draws  its  name  from  the  manner  in  which  the  search 
direction,  s,  is  determined.  If  the  current  point  xk  is  on  the  inside  of  the  feasible  region  the 
search  direction  is  the  direction  of  the  gradient  of  the  objective  function.  If  the  current  point 
is  on  the  border  of  the  feasible  region,  then  the  search  direction  is  the  projection  of  the 
gradient  on  the  binding  constraint(s).  This  results  in  the  search  for  the  optimum  point 
moving  along  the  binding  constraints  in  the  second  case.  The  second  case  is  illustrated  in 

Figure  1.  Active  constraints  are  those  for  which  the  equality  holds  -  i.e.  gi(xk)  =  cit  and 
h  (x  J  =  b  .  Binding  constraints  are  those  active  constraints  which  restrict  the  direction  of 

the  optimization. 
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Figure  1:  Projection  of  the  gradient  of  the  objective  function  on  the  binding  constraint 

The  projection  matrix  is  calculated  using 

P  =  [  I  -  AyAA^A],  (C.2) 

where  the  A  matrix  consists  of  the  coefficients  of  the  binding  constraints.  Each  binding 
constraint  is  a  row  in  A.  In  order  to  determine  which  of  the  active  constraints  are  binding 
the  La  Grange  multipliers,  X,  can  be  used.  The  multipliers  are  calculated  according  to 

X  =  A‘(AaA‘)'lvf(xk)’  (C.3) 

where  Aa  is  the  matrix  of  all  active  constraints  with  each  active  constraint  occupying  a  row 
of  Aa.  The  positive  elements  of  X  indicate  that  the  corresponding  constraints  are  binding. 
The  search  direction  is  then  calculated  as 


s  =  PVf.  (C.4) 

The  incorporation  of  the  nonlinear  constraints  causes  some  difficulty.  In  order  to 
incorporate  the  nonlinear  constraints  they  are  linearized  to  be  the  gradient  of  the  constraint 
at  the  point  of  interest.  The  projection  matrix  is  calculated  using  these  linearized  constraints. 
This,  however,  creates  problems. 
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Let  y  =  Xk  +  as  where  x^  is  a  feasible  point,  s  is  a  direction  which  violates  none  of  the 
linear  or  linearized  constraints,  and  a  is  the  step  length.  Since  the  linearized  constraints  are 
only  accurate  at  the  point  the  point  y  will  violate  the  nonlinear  constraints.  This  is 
demonstrated  in  Figure  2. 


Figure  2:  Projection  of  the  gradient  of  the  objective  function  on  a  linearized  binding 
constraint 

In  order  to  return  to  the  feasible  region  the  the  following  iterate  process  is  used. 

yj+i  =  yj  *  AHAA'j'kgbCyj)  -  cb]  c.5 

where 

Sb(yj)  is  the  set  of  the  all  of  the  binding  constraints  evaluated  at  yj  and  cb  is  the  set  of 
corresponding  boundary  conditions.  This  process  will  never  return  the  point  exactly  to  the 
feasible  region,  but  it  can  come  to  within  8  of  the  feasible  region.  Let  w(y)  refer  to  the 
results  of  the  process  for  the  infeasible  point  y. 

Outline  of  the  Algorithm 


The  basic  idea  of  the  optimization  algorithm  is  to  move  from  the  current  feasible  point  to  a 
better  point.  This  process  of  moving  from  point  to  point  in  the  direction  can  be  broken 
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down  into  two  steps:  1)  calculating  the  search  direction  and  2)  determining  how  far  to  move 
in  the  search  direction  to  the  next  point. 

STEP  0 

start  at  an  initial  feasible  point  xo,  and  choose  a  maximum  step  for  the  nonlinear  case.  The 
nonlinear  step  size  is  denoted  by  «n 

STEP  1:  Calculating  the  search  direction  s 

The  first  step  can  be  broken  down  into  several  parts: 

1.  Calculate  Vf(xk),  the  gradient  of  f(x^). 

2.  Determine  the  active  constraints,  and  find  the  binding  ones  from 
X  =  A'CAXl'Vf^)  f  X,  >  0 

(Note  that  most  algorithms  don't  differentiate  between  binding  and  nonbinding 
active  constraints  until  s  =  0.  This  is  done  for  computational  reasons.) 

3.  Use  the  binding  constraints  to  create  A.  The  projection  matrix  can  then  be 
calculated  according  to 

P  =  [  I  -  AKAA^A] 

(Note  that  if  none  of  the  constraints  are  binding  P  =  I,  the  identity  matrix.) 

4.  The  search  direction  is  then  calculated  from 
s  =  PVf 

If  the  absolute  value  of  the  search  direction  is  zero  then  the  constrained  optimum 
has  been  reached.  For  practical  purposes  the  optimum  is  considered  to  have  been 
reached  if  the  absolute  value  of  s  is  sufficiently  small. 
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5.  if  |s|  <  e  then  stop. 

Step  2:  determining  how  far  to  move  in  the  search  direction  to  the  next  point 

The  second  step  is  broken  down  into  two  cases.  The  first  case  is  when  the  nonlinear 
constraints  are  nonbinding.  The  second  is  when  the  nonlinear  constraints  are  binding. 

Case  1 :  The  nonlinear  constraints  are  not  binding 

1.  First  the  maximum  step  size  is  determined. 

Max  {  <Xj:  u  =  x  +  otjS  is  feasible  } 

2.  The  best  point  along  this  line  is  determined.  This  line  search  is  done  using  a 
secant  search.  (This  is  a  one  dimensional  search.)  In  the  code  for  this  thesis  the 
algorithm  stops  after  ni  consecutive  line  searches. 

Max  {  a2:  f(xk  +  a2s),  0  <  a2  <  otj  } 

3.  The  current  point  is  then  updated  to  be  this  point. 
xk+i  =  xk  +  a2s 

4.  return  to  Step  1 

Three  examples  of  Step  1  and  Case  1  of  Step  2  are  shown  in  Figures  3,4,5. 
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Figure  3:  The  search  direction  is  the  direction  of  the  objective  function  gradient.  The 
optimum  point  in  this  direction,  x^+i^  is  between  the  bound  and  the  current  point,  xj<. 


Figure  4:  The  search  direction  is  the  direction  of  the  objective  function  gradient.  There  are 
two  active  constraints  at  the  point  x^  but  neither  is  binding.  The  optimum  point  in  this 
direction,  xfc+i,  is  on  the  bound. 
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Figure  5:  The  search  direction  is  the  direction  of  the  projection  of  objective  function 
gradient  on  the  binding  constraint.  The  optimum  point  in  this  direction,  x^+i^  is  on  the 

bound. 

Case  2:  at  least  one  of  the  nonlinear  constraints  is  binding 

If  the  any  of  the  nonlinear  constraints  are  binding  then  a  new  procedure  must  be  followed. 
This  procedure  consists  of  taking  steps  of  size  an  along  the  nonlinear  constraints  as  long  as 
the  step  results  in  an  improved  point.  The  step  size  an  is  halved  if  the  new  point  does  not 
result  in  a  better  point. 

1 .  Find  a  point  where  the  projection  results  in  a  feasible  solution.  For  the  code 
written  for  this  thesis,  this  is  done  by  cutting  the  step  length  «n  in  half  until  a  point 
is  found  where  the  projection  back  on  to  the  constraints  results  in  a  feasible  point 
(see  Figure  6). 

Min  {p:  w(y  =  xjc  +  (1/2)pOjiS)  is  feasible} 
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Vf(xk) 


*2 


Xl 

Figure  6:  Halving  step  size  to  Find  a  feasible  projection  on  constraints 

2.  Next  determine  if  this  new  point  resulted  in  an  increase  in  the  objective  function. 

a)  If  it  resulted  in  an  increase,  then  update  the  current  point  to  be  this  new 
point  and  the  step  size  remains  the  same  (see  Figure  7).  In  order  to  speed 
the  optimization  process,  the  step  size  is  doubled  if  seven  consecutive  steps 
at  the  same  step  size  have  resulted  in  improved  points  and  doubling  the  step 
size  will  result  in  a  step  size  less  than  the  maximum  step  size  for  the 
nonlinear  case  set  in  Step  0. 

If  f[w(y)]  £  f(xk)  then  a£+1  =  a£  and  xk+!  =w(y). 
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Figure  7:  Current  point  is  updated  to  be  w(y),  the  point  on  the  nonlinear  constraint  when  y 
is  a  better  point  than  x k 

b)  If  the  new  point  is  an  inferior  point,  the  step  size  is  cut  in  half  and  the 
current  point  remains  the  same. 


If  f[w(y)J  <  f(xk) ,  then  a£+1  =  l/2a£  and  xk+J  =  xk. 


3.  If  otn  -  e2  or  an  <  £3,  f(xk+1)  -  f(xk)  >  e4  and  the  number  of  cycles  is  greater 
than  nj’  t^en  stoP’  otherwise  return  to  step  1. 
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APPENDIX  D 

THE  SEQUENTIAL  OPTIMIZATION  PROGRAM 


#include  <stdio.h> 

#include  <math.h> 

y******************************************************************y 

/*  global  variables  */ 

y******************************************************************y 

double  v[50][50],kl,bc[20],con_mat[50][50],cur[20],aaa[50][50],sigma; 
int  pointer,ncon,flag,novar,ncoeff,pl,flag2,flag3,flag9; 

double  c[30],g[20],s[20],t,aa[50][50],sig2[20],ave[20],hv[10]; 

double  var2[50][50],ave2[20],d; 

y******************************************************************y 
main(argc,argv) 
int  argc; 
char  *argv[]; 

{ 

y* ******************************************************************y 

/*  external  variables  in  main  */ 

y*******************************************************************y 

extern  double  v  [  50]  [ 50]  ,k  1  ,bc[20]  ,con_mat[50]  [50]  ,cur[20]  ,aaa[50]  [50]  ,sigma; 

extern  int  pointer,ncon,fiag,novar,ncoeff,flag2; 

extern  double  c[30],g[20],s[20],t,sig2[20],ave[20],hv[10]; 

y* ********************************************************** **********y 

/*  functions  in  main  */ 

y* ********************************************************** **********y 

double  abss(),agrad(),bk(),bound(),create(),conck(),conse(),cpm(),grad(); 
double  lines(),nim(),dconck(),dconse(); 

double  modelg(),mvm(),nran(),reg(),trans(),vran(),w(),wbound(); 

y* ******************************** ***********************************y 

/*  declaration  of  variables  in  main  */ 

y********************************************************************y 

PILE  *sp,*sp3,*sp9; 

double  b[50] [50] ,a[50]  [50] ,proj [20] , slope, f  1  ,f3  ,alphan; 
double  y,at,alpha_old; 
double  xtr[50][50]; 

int  check,i,j,l,nodat,k,r,p,timeold, repeat; 
double  lv[10],ad,time,  fI20]junk; 
double  curm[20],ym,dck; 
static  double  st[32]  =  { 

[0,6.3,2.9,2.3,2.1,2.0,1.94,1.9,1.86,1.83,1.812,1.8,1.78,1.77,1.76,1.75, 

1.75,1.74,1.73,1.73,1.725,1.72,1.72,1.71,1.71,1.71,1.7,1.7,1.7,1.7,1.7,1.64], 

); 

static  double  chi[17]={0.0, 2.71, 4.61, 6.25,7.78, 9.24,10.64,12.02,13.36,14.68, 
15.99,17.28,18.55,19.81,21.06,22.31,23.54); 
double  alpha,pm[50][50],sub(); 
double  ab,am,tck,yhat; 

/*  random  variable  seed  is  set  */ 

srand(10); 

y********** ************************************************** **********y 
/*  beginning  of  algorithm  */ 

y**********************************************************************y 
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for(p=l  ;p<=25:++p) 
r 

i 


if(p==l) 

{ 

/********************************************+***********************/ 


I*  beginning  of  initialization  routine  */ 

/*  this  step  is  only  executed  before  the  first  cycle  */ 

f*  it  asks  the  user  to  set  the  parameters  such  as  the  scaling  */ 

/*  constants  and  the  initial  reference  value  of  y,  Kl=  .l*yref  */ 


/*********4^************************%***%****%*****%*****************y 


printf("how  many  variables?\n"); 
scanf("%d",&novar); 
printf("what  is  y  refNn"); 
scanf("%f',&y); 
printf("what  is  #  coefficientsNn”); 
scanf("%d”,&ncoeff); 
t=pow(y*y,0.5)*0. 1 ; 
for(j = 1  ;j  <=nov  ar;  ++j ) 

{ 

printf("what  are  high  and  low  values  for  variable  %d?\n"  j); 

scanf("%f  %f ’,&hv[j],&lv[j]); 

hv[j]=hv[j]-lv[j]; 

} 

sp3=fopen("numdat","r"); 

fscanf(sp3,"%d",&nodat); 

fclose(sp3); 


} 

^*  ******************+*******>*■************+****************>»<**»******  *^/ 
/*  advice  is  given  and  the  results  of  the  experiment  are  collected  */ 

/****************************************************.****************/ 

if(p>=2) 

{ 


vran(cur,curm); 
printfC’run  expt  at\n"); 
for(j=l  ;j<=novar,++j) 

{ 

printf( "variable  %d  =  %f\n",j,cur[j]); 

} 


/*  input  new  data  and  update  numdat  and  exper.dat  files  */ 

printf("what  is  yTVn”); 

scanf("%f',&y); 

ym=y/*+nran()*0.005;*/; 

printf("y  measured  =  %f\n",ym); 

sp3=fopen("numdat","r"); 

fscanf(sp3,"%d",&nodat); 

fclose(sp3); 
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nodat=nodat+l; 
sp3=fopen("numdat"  ,"w"); 
fprintf(sp3,"%d\n",nodat); 
fclose(sp3); 
t=0.1*pow(y*y,0.5); 

/*  append  new  data  to  exper.dat  file  */ 

printfO'enter  the  data  actually  runW); 
for(k=  1  ;k<=novar;++k) 

{ 

scanf("%f',&cur[k]); 
printfC’got  the  data\n"); 

) 

sp=fopen("exper.daf,"a"); 
fprintf(sp,"%f  ",ym); 

for(k=  1  ;k<=novar;++k) 

{ 

fprintf(sp,"%f  ”,cur[k]); 

} 

fprintf(sp,"\n"); 

fclose(sp); 

/******************************************************************/ 

/*  book  keeping  */ 

I ** ** *** ** * ** * ** *****  *  ****** ** ****** ** ** ** *****  ******************* j 

yhat=modelg(cur); 

tck=conck(cur); 

bk(yhat,p,tck,y); 

} 


/*********************************************************************/ 
/*  call  reg  to  update  the  coefficients  */ 

I*********************************************************************/ 

reg(nodatjiovai\hv,v,&sigma,g,c); 

I*********************  ***********************************************/ 

/*  calculate  the  kl  value  (the  appropriate  t)  */ 

/*********************************************************************/ 
d=2.0*chi[novar]; 
r=nodat-ncoeff; 
if  (r>=31) 

{ 

r=31; 

) 

kl=st[r]; 


j=0; 

for(l=  1  ;l<=novar,++l) 

{ 

cur[l]=g[l]; 

) 
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time=0.0; 

check=0; 

repeat=0: 

alpha_old=0.0; 

yT*  ********************************%*************%***%*/ 

/*  start  of  optimization  */ 

/**********************»******************************/ 
at  =1.0; 
while  (j<=0.0) 

{ 

/*  first  the  search  direction,  s.  is  determined  */ 

grad(f); 

cpm(pm,f,time); 

mvm(s,pm,f,novar,novar); 

/*  if  Isl  is  small  then  the  optimization  is  completed  */ 

if(abss(s,novar)<=0.00 1 ) 

( 

if(time==0.0) 

{ 

j=i; 

printf("help\n"); 
goto  end; 

} 

j=i; 

tck=conck(cur); 
printf("  tck  =  %f\n",tck); 
dck=dconck(cur); 
printf("  dck  =  %f\n",dck); 

} 

else 

{ 

/*  case  where  neither  of  the  nonlinear  constraints  is  binding  */ 

if  (flag==0  &&  flag3==0) 

{ 

/*  finds  the  boundary  that  will  first  be  encountered  */ 

/*  by  travelling  in  the  direction  of  the  s  vector  */ 

ab=bound(ncon,novar,s,cur,con_mat,bc, alpha); 

am  =conse(xtr,g,novar,a,b,st,r,ncoeff,v,sigma,ab,s,t,cur); 

ad  =dconse(xtr,g,novar,a,b,st,r,ncoeff,v,sigma,ab,s,t,cur); 

alpha=am; 

if(am>=ab) 

{ 

alpha=ab; 

} 

if(ad<alpha) 

{ 
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alpha=ad; 


/*  if  slope  is  positive  at  the  bound  then  next  point  */ 

/*  is  point  at  the  bound  do  line  search  in  direction  of  */ 
/*  the  s  vector  */ 

slope=agrad(alpha); 

if(slope>0.0) 

{ 

for  (i=l;i<=novar;++i) 

{ 

cur[i]=cur[i]+alpha*s[i]; 

} 

) 

else 

/*  if  slope  is  negative  at  bound  then  a  secant  line  search  is  */ 

/*  performed  by  the  subroutine  lines  to  find  nest  point  */ 

{ 

junk=lines(0.0,alpha, slope); 
if(junk  <=  0.000000001) 

{ 


j=i; 

tck=conck(cur); 
printf("  tck  =  %f\n",tck); 
dck=dconck(cur); 
printf("  dck  =  %f\n",dck); 

} 

for  (i=l;i<=novar;++i) 

{ 

cur[i]=cur[i]+junk*s[ij; 

} 

/*  check  #  of  line  searches,  quit  if  it  is  to  large  */ 

if(time=timeold+ 1 ) 

{ 

repeat=repeat+l; 

} 

else 

{ 

repeat=0; 

} 

if(repeat>=50) 

{ 

j=i; 

tck=conck(cur); 

repeat=0; 

printfO’not  repeating  another  line  search"); 

} 

timeold=time; 

) 
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else 


/*  the  nonlinear  constraints  are  binding  */ 


alpha=wbound(at,proj); 
if  (alpha>=at)  alpha=at; 

/*  counting  #  of  consecutive  steps  at  a  given  size  */ 

if(alpha==alpha_old) 

{ 

check=check+l; 
if  (check>=8)  check=0; 
if(check>=7  &&  alpha<=0.5) 

{ 

alpha=2.0*alpha; 

at=alpha; 

} 

} 

else 

{ 

alpha_old=alpha; 

check=0; 

} 

/*  evaluating  the  model  at  old  and  new  point  */ 

fl=modelg(cur); 

f3=modelg(proj); 

/*  printf("  f  1  =  %{,  f3=%f\n”,fl  ,f3);*/ 

/*  update  current  point  to  be  this  new  point  if  new  point  is  better  */ 

if(f3>=fl) 

{ 

for  (i=l;i<=novar;++i) 

{ 

cur[t]=proj[i]; 

} 

junk=B-fl; 

/*  check  for  stopping  conditions  */ 

if(at<=.001  &&  time>=1500  &&  junk<=0.0000001) 

( 

printf("we're  taking  the  coward's  route  home”); 

tck=conck(cur); 
printf("  tck  =  %f\n",tck); 
dck=dconck(cur); 
printfC  dck  =  %f\n”,dck); 
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j=i; 

} 


/*  if  new  point  is  better  then  halve  step  size 
else 
{ 

at=0.5*at; 

} 

/*  check  for  stopping  conditions 

if  (alpha  <=  0.0000001) 

{ 

j=i; 

tck=conck(cur); 
printf("  tck  =  %f\n",tck); 
dck=dconck(cur); 
printf("  dck  =  %f\n'',dck); 

} 

} 

} 

time  =dme+l; 


*/ 


} 

j*  ****************************************************************^ 

/*  end  of  optimizer  */ 

j *  **************************************************************  *  j 

} 

end:  /*  this  is  it  */; 

} 

/** *************************************** **********************/ 


/* 


end  of  main  program  and  start  of  subroutines 
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*/ 

*/ 


/***************************************************************/ 

/*  subroutine  to  do  matrix  multiplication 

/*  a=xl*x2  xl  is  nxp  and  x2  is  pxm 

/*  ***************************************************** *********/ 
double  mm(a,xl,x2,n,p,m) 
int  n,m,p; 

double  a[50][50],xl[50][50],x2[50][50]; 

{ 

int  i,j,l; 

for  (i=l;i<=n;++i) 

{ 

for(j=l;j<=m;++j) 

{ 

a[i][j]=0-0; 

for(l=l;l<=p;++l) 

{ 

a[i][j]=a[i][j]+xl[i][l]*x2[l][j]; 

} 

) 

} 

} 

y* ******************************************** *********************/ 
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/*  subroutine  to  create  predictor  variables  from  control  variable  */ 

/******************************************************************/ 

double  create(x,xr,nodat) 

double  xr[20],x[50][50]; 

int  nodat; 

l 

extern  int  flag2,flag9; 
int  j,l,  counter, 
double  temp[20]; 
for(j=l;j<=3;++j) 

{ 

temp[j]=x[nodat]|j]; 

} 

/*  check  to  see  if  conversion  will  cause  math  error,  if  so  flag2  =  0  */ 

if(temp[l]<0.0  II  temp[3]<=0.0) 

{ 

flag2=0; 
goto  qu; 

} 

flag9=l; 

x[nodat][l]=1.0; 

x[nodat][2]=temp[l]; 

x[nodat][3]=temp[2]; 

x[nodat][4]=temp[3]; 

x[nodat][5]=temp[  1  ]*temp[  1  ]; 

x[nodat][6]=temp[2]*temp[2]; 

x[nodat][7]=temp[3]*temp[3]; 

x[nodat][8]=temp[l]*temp[2]; 

x[nodat][9]=temp[l]*temp[3]; 

x[nodat][10]=temp[2]*temp[3]; 

qu:; 

}* 

j *  ♦ilc***************************************************************^ 

/*  subroutine  to  transpose  x-vector  */ 

/*  xt  is  the  transpose  of  x,  x  is  lxm  vector  */ 

/♦★♦♦★fr*************************************************************/ 

double  trans(x,xt,p) 

double  x[50][50],xt[50][50]; 
int  p; 

{ 

int  m; 

for(m=l  ;m<=p;++m) 

{ 

xt[m][l]=x[13[m]; 

) 

} 

j*  **************************************★************************/ 

/*  subroutine  to  multiply  matrix  and  vector  */ 

/*  a  =  h*vec  hisnxp,  vecispxl  */ 

/****************************************************************/ 
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double  mvm(a,h,vec,n,p) 

double  a[20],h[50][50],vec[20]; 
int  n,p; 

{ 


int  j.i; 

for  (i=l;i<=n;++i) 

{ 

a[i]=0.0; 

for  (j=l;j<=p;++j) 

{ 

a[i]=a[i]+h[i][j]*vec[j]; 

} 


} 

/* 

/* 

/* 

/* 


} 


* ************************************************************* 


gradient  for  prediction  error  constraint 
is  numerically  determined 

************************************************************** 


double  gradh(f.vec) 
double  f[20],vec[20]; 


*/ 

*/ 

*/ 

*/ 


extern  novar,ncoeff; 
int  l,i; 

double  conck(),xp[20],p[20]; 
for(l=0;  l<=novar,  ++1) 

{ 

xp[l]=vec[l]+.00000001; 

for(i=l;i<=l*l;++i) 

{ 

p[i]=vec[i]; 

} 

if(l>=l) 

{ 

p[l]=xp[l]; 

} 

for(i=l+ 1  ;i<=novar,+-fi) 

{ 

p[i]=vec[i]; 

} 

f!l]=conck(p); 
if  (1>=1) 

{ 

f[l]=(fll]-f[0])/.00000001; 
f*  printf("%f\n",f[l]);*/ 

} 

} 

} 

y*  ********************************************************** *****/ 
/*  gradient  for  distance  constraint  */ 

/*  is  numerically  determined  */ 

I*  ******************************************** ******************  *j 

double  dgradh(f,vec) 
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double  f[20],vec[20]; 


extern  novar,ncoeff; 
int  l,i; 

double  dconck(),xp[20],p[20]; 
for(l=0;l<=novar,++l) 

{ 

xp[l]=vec[l]+.00000001; 

for(i=l;i<=l-l;++i) 

{ 

p[i]=vec[i]; 

} 

if(l>=l) 

{ 

p[l]=xp[l]; 

} 

for(i=l+l;i<=novar;++i) 

{ 

p[i]=vec[i]; 

} 

f!l]=dconck(p); 
if  (1>=  1 ) 

{ 

f[l]=(fll]-f[0])/.00000001; 

/*  printf(”%f\n",f[l]);*/ 

} 

} 

} 

I*  ************************************************************** 

/*  gradient  for  objective  function  */ 

/*  is  numerically  determined  */ 

/**************************************************+*************/ 
double  grad(0 
double  f[20]; 

{ 

extern  double  cur[20]; 
int  l,i; 

double  modelg(),xp[20],p[20]; 
for(l=0;l<=novar,++l) 

{ 

xp[l]=cur[l]+.0000000 1 ; 
for(i=  1  ;i<=l- 1  ;++i) 

{ 

p[i]=cur[i]; 

} 

if(l>=l) 

{ 

p[l]=xp[l]; 

} 

for(i=l+l  ;i<=novar,++i) 

{ 

p[i]=cur[i]; 
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} 

fIl]=modelg(p); 
if  (1>=1) 

{ 

f[l]=(f[l]-f[0])/. 00000001 ; 

/*  printf("%f\n”,f[l]);*/ 

} 

} 

} 

I*  *********************************************************  * j 

/*  subroutine  to  create  p  matrix  */ 

/*  the  P  matrix  is  a  projection  matrix  */ 

j *  *********************************************************  *i 

double  cpm(pm,f,time) 

double  pm[50][50],fI20],time; 

{ 

extern  double  cur[20],  con_mat[50][50],aaa[50][50],bc[20],g[20]; 
extern  int  pointer,novar,coeff,flag,flag3; 
extern  double  aa[50][50],kl,t,d; 

FILE  *sp2; 

int  i,l,j,m,junk,count,info; 

double  aat[50]  [50]  ,a2[20]  ,a[50]  [50]  ,b[50]  [50]  ,fh[20] ; 

double  dgedi_(),dgefa_(),mvm(),ipvt[50],work[50],det[3],u[20]; 

double  conck(),dconck(),gradh(),dgradh(); 

sp2=fopen("constr","r"); 

fscanf(sp2,"%d",&ncon); 

/*  initializing  the  pm  matrix  by  setting  it  equal  to  0  */ 
for(i=l  ;i<=novar,++i) 

{ 

for(l=  1  ;l<=novar,++l) 

{ 

pm[l][i]=0.0; 

} 

} 


/*  read  in  user  constraints  from  constr 

for  (i=l;i<=ncon;++i) 

{ 

for(l=  1  ;l<=novar,++l) 

{ 

fscanf(sp2,"%f',&con_mat[i][l]); 

} 

fscanf(sp2,"%f ',&bc[i] ); 


} 

fclose(sp2); 

/*  set  up  aa  matrix  with  active  constraints 


*/ 


*/ 


mvm(a2,con_mat,  cur,  neon,  novar); 

pointer=0.0; 

for  (i=l;i<=ncon;++i) 
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if(bc[i]-a2[i]>=-.01) 

{ 

pointer=pointer+l : 
for(m=l  ;m<=novar:++m) 

{ 

aa[pointer][m]=con_mat[i][m]; 


} 


} 


/*  see  if  nonlinear  prediction  constraint  is  active  */ 

if(conck(cur)>=t-.01 ) 

{ 

if(time==0) 

{ 

t=conck(cur); 
printfC'changing  t\n"); 

} 

pointer=pointer+ 1 ; 
gradh(fh,cur); 
for(i=  1  ;i<=novar,++i) 

{ 

aa[pointer][i]=0.0  -fh[i]; 

} 

flag=l; 

} 

else 

{ 

flag=0; 

} 

/*  (flag  and  flag3  are  global  variables  that  indicate  whether  or  not  the*/ 

/*  prediction  error  constraint  and  the  distance  constraint  respectively  are*/ 
/*  binding)  */ 

/*  check  to  see  if  the  nonlinear  distance  constraint  is  active  */ 

if(dconck(cur)>=d-.0 1 ) 

{ 

if(time==0) 

( 

d=dconck(cur); 
printfC'changing  dNn"); 

} 

pointer=pointer+l; 

/*  if  the  constraint  is  active  linearize  the  nonlinear  distance  */ 

/*  bound  by  using  its  gradient  */ 

dgradh(fh,cur); 
for(i=  1  ;i<=novar,  ++i) 

{ 
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aa[poimer][i]=0.0  -fh[i]; 

) 

flag3-l; 

} 

else 

{ 

flag3=0; 

} 

for(j=l  ;j<=pointer;++j) 

{ 

for(i=l;i<=novar,++i) 

{ 

aat[i][j]=aa[j][i]; 

} 

} 

/*  set  up  the  identity  matrix  */ 

if(pointer  ==  0) 

{ 

for  (i=l;i<=novar;++i) 

{ 

pm[i][i]=1.0; 

} 

} 

else 

( 

for  (i=l;i<=novar;++i) 

{ 

pm[i][i]=1.0; 

} 

mm(a,aa,aat,pointer,novar,pointer); 

/*  invert  matrix  */ 
for(i=  1  ;i<=pointer;++i) 

{ 

for(m=  1  ;m<=pointer;++m) 

{ 

a[i*l][m-l]=a[i][m]; 

} 

} 

dgefa_(a,&pointer,&pointer,ipvt,&info); 
dgedi_(a,&pointer,&pointer,ipvt,det,work,01); 
for(i=pointer- 1  ;i>=0; — i) 

{ 

for(m=pointer- 1  ;m>=0;--m) 

{ 

a[i+l][m+l]=a[i][m]; 

} 

} 

mm(b,a,aa, pointer, pointer, novar); 

/*  check  LaGrange  multipliers  to  see  which  active  constraints  are  binding  */ 
mvm(u,b,f,pointer,novar); 
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for(i=2;i<=pointer+ 1  ;++i) 

( 

count=0; 

for(l=i- 1  ;1>=  1 ; — 1) 

{ 

if(u[l]>0.0) 

{ 

count=count+l; 

if(l==pointer-flag3  &&  flag==l) 

{ 

flag=0; 

printf("changed  flag\nM); 

) 

if(l==pointer  &&  flag3=— 1 ) 

{ 

flag3=0; 

printf("changed  flag3\n"); 

) 


} 

} 

if(count>0) 

{ 

for(j = 1 ; j  <=novar, ++j ) 

{ 

aa[i-count]  [j]  =aa[i]  [j] ; 

} 

} 

} 

/*  pointer  is  global  variable  and  is  #  of  binding  constraints  */ 

pointer=pointer-count; 

/*  construct  projection  matrix  */ 

for(j=l  ;j<=pointer;++j) 

{ 

for(i=l  ;i<=novar,++i) 

( 

aat[i][j]=aa[j][i]; 

} 

} 

mm(a,aa,aat,pointer,novar,  pointer); 
f*  invert  matrix  */ 
for(i=  1  ;i<=pointer;++i) 

{ 

for(m=l  ;m<=pointer;++m) 

{ 

a[i-l][m-l]=a[i][m]; 

} 

) 

dgefa_(a,&pointer,&pointer,ipvt,&info); 
dgedi_(a, &po  inter, &pointer,ipvt,det,work,01); 
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for(i=pointer- 1  ;i>=0;--i) 

{ 

for(m=pcmter-l;m>=0;-m) 

{ 

a[i+l][m+l]=a[i][m]: 

} 

) 

mm(aaa,aat,a,novar,pointer,pointer); 

mm(a,aaa,aa,novar,pointer,novar); 

for(i=l;i<=novar;++i) 

{ 

for(l=l;l<=novar;++l) 

{ 

pm[i][l]=pm[i][l]-a[i][l]: 

} 

} 


j *  ***************************************************  */ 

/*  alpha  grad  numerically  calculates  df/d(alpha)  */ 

j *  *  *  **  **  ***  *  *  *  ***  **  **  *  *  ***  *  *****  ***  *******  **  **  J(c  *  *  *  *  *  * 

double  agrad(alpha) 
double  alpha; 


{ 


int  l,i; 

double  modelg(),grad[3],result,xp[20],junk; 
for(l=0;l<=l;++l) 

{ 

w(.0000 1  *l+alpha,xp); 
grad[l]=modelg(xp); 

} 

result=(grad[  1  ]-grad[0]  )/.0000 1 ; 

/*  printf("%f  ".result);*/ 
return  (result); 

) 

jUf  **+*****************************************+****+*  */ 

/*  bound  finds  the  stepsize  in  the  s  direction  */ 

/*  to  the  nearest  (linear)  user  constraint  */ 

j*  ***************************************************  *i 

double  bound(ncon,novar,s,g,h, be, alpha) 
double  h[50][50],s[20],g[20],bc[20], alpha; 
int  novar,ncon; 

{ 

double  mvm(),atx[20],ats[20],amax[20],fabs(); 
int  i; 

mvm(atx,h,g,ncon,novar); 

mvm(ats,h,s,ncon,novar); 
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/*  printfC  %i  %{  %fNn",s[l],s[2],s[3]);*/ 
for(i=  1  ;i<=ncon;++i) 

{ 


if(ats[i]==0.0  II  fabs(ats[i])  <=  l*pow(10.0,-7.0)) 

{ 

am  ax[i]=  1000.0; 


} 


} 

else  if  (bc[i]/ats[i]-atx[i]/ats[i]>0.0) 

{ 

amax[i]=(bc[i]-atx[i])/ats[i]; 

) 

else 

{ 

amax[i]= 1000.0; 

} 


/*  find  alpha  */ 
alpha=  1000.0; 
for(i=l  ;i<=ncon;++i) 

{ 

if(amax[i]<=alpha) 

{ 

alpha=amax[i]; 

} 

} 

retum(alpha); 


} 

/ft********************************************************/ 

/*  conck  evaluates  1/2  the  width  of  the  prediction  */ 

/*  confidence  interval  at  x  */ 

/*************************#*+*****************************/ 
double  conck(x) 
double  x[20]; 


{ 

extern  double  g[20],v[50][50], sigma, kl,hv[10],var2[50][50],ave2[20]; 

extern  int  novar.ncoeff.nodat; 

inti; 

double  xt[50][50],xtr[50][50],a[50][50],b[50][50],create(),trans(),mm(),standard0; 
double  w,c[50][50],f; 
for(i=l  ;i<=novar,++i) 

{ 

xtr[l][i]=x[i]; 

) 

w=0.0; 


for(i= 1  ;i<=novar, -H-i) 

{ 

w=w+pow((x[ij-g[i])/hv[i],2.0); 


} 

w=w/novar/2.0; 

if(w>=10.0) 

{ 
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w=10.0; 

} 

w=exp(w); 

/*  test  for  distance  */ 
for(i=l;i<=novar,++i) 

{ 

c[i][l]=x[i]-ave2[i]; 

} 

mm(a,var2,c,novar,novar,  1 ); 
for(i=l;i<=novar,++i) 

{ 

c[l][i]=c[i][l]; 

} 

mm(b,c,a,  1  ,novar,  1 ); 
f=b[  1][1  ]; 
create(xtr,g,l); 
standard(xtr); 
trans(xtr,xt,ncoeff); 
mm(b,xtr,v,  1  ,ncoeff,ncoeff); 
mm(a,b,xt,  1  .ncoeff,  1 ); 

a[l][l]=pow((w+a[l][l])*sigma,.5)*kl; 

/*a[l][l]=pow(a[l][l],.5)*kl;*/ 

retum(a[l][l]); 

} 

/**************************************************£*« *************/ 
/*  dconck  evaluates  the  distance  at  x  */ 

/******]^**** ******************************************** **********/ 
double  dconck(x) 
double  x[20]; 


{ 

extern  double  g[20],v[50] [50],sigma,k  1  ,hv[l 0],var2[50] [50] ,ave2[20]; 
extern  int  novar,ncoeff,nodat; 
int  i; 

double  xt[50][50],xtr[50][50])a[50][50],b[50][50],create(),trans(),mm(),standard(); 
double  w,c[50][50],f; 

/*  test  for  distance  */ 
for(i= 1  ;i<=novar,++i) 

{ 

c[i]  [  1  ]=x[i]-ave2[i]; 

} 

mm(a,var2,c,novar,novar,  1 ); 
for(i=l  ;i<=novar,++i) 

{ 

c[l][i]=c[il(l]; 

} 

mm(b,c,a,l,novar,l); 
f— b[l][l]; 
retum(f); 

} 

^*  ******************************************************  */ 
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/*  conse  finds  the  maximum  step  size  based  on  the  */ 

/*  prediction  error  constraint  */ 

/*  ****************************************************** 
double  conse(xtr,g,novar,a,b,st,r.ncoeff,v,sigma,alpha,s,t,cur) 
double  alpha,xtr[50] [50], g[20],a[50][50],b[50][50],st[32], sigma; 
double  v[50][50],s[20],t,cur[20]; 
int  novar,r,ncoeff; 

{ 

extern  double  hv[10]; 
int  kj,i; 

double  alphau,alphal,xl[50][50],x2[50][50],xt[50][50],w; 

double  create(),standard(); 

for(k=l;k<=novar;++k) 

{ 

xl[l][k]=cur[k]+alpha*s[k]; 

} 

for(k=  1  ;k<=novar;++k) 

{ 

x2[l][k]=cur[k]: 

} 

j=0; 

alphau  =  alpha; 
alphal=0.0; 

/*  if  the  calculated  variance  of  y-predicted  is  not  less  than  maximum  */ 

/*  allowable  value  iterate  to  find  a  point  equal  to  the  maximum  */ 

while(j<=0) 

{ 


alpha=0.5  *(alphau+alphal); 
for(k=l  ;k<=novar;++k) 

{ 

xtr[  1  ]  [k]=cur[k]+s[k]  *alpha; 

) 

w=0.0; 

for(i=  1  ;i<=novar,++i) 

{ 

w=w+pow((xtr[  1  ]  [i]-g[i])/hv[i],2.0); 

) 

w=w/2.0/novar; 

if(w>=10.0) 

{ 

w=10.0; 

} 

w=exp(w); 
create(xtr,g,l); 
standard(xtr); 
trans(xtr,xt,ncoeff); 
mm(b,xtr,v,  1  ,ncoeff,ncoeff); 
mm(a,b,xt,  1  .ncoeff,  1 ); 
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a[l][l]=pow((w+a[l][l])*sigma,.5)*st[r]; 

if(a[l][l]>=t+.0001) 

{ 

for(k=l;k<=novar;++k) 


alphau=alpha; 

) 


} 

if(a[  1  ]  [  1  ]  <=t- .000 1 ) 

{ 


for(k= 1 ;  k<=novar; ++k) 

{ 


alphal=alpha; 


} 

} 

/*  condition  for  exiting  the  while  loop  */ 
if(a[l][l]<=t+.0001  &&  a[l][l]>=t-.0001) 
{ 

for(k=  1  ;k<=novar;++k) 

{ 

xtr[l][k]=xtr[l][k]; 

} 

j=i; 


if(alphau-alphal<=.0000 1 )  j= 1 ; 

} 

retum(alpha); 

} 

/*  ****************************************************** 

/*  dconse  finds  the  maximum  step  size  based  on  the 

/*  distance  constraint 

j *  ****************************************************** 


*/ 

*/ 

*/ 

*/ 


double  dconse(xtr,g,novar,a,b, st, r,ncoeff,v, sigma, alpha, s,t,cur) 
double  alpha,xtr[50] [50] ,g[20] ,a[50] [50] ,b[50] [50] ,st[32] , sigma; 
double  v[50][50],s[20],t,cur[20]; 
int  novar,r,ncoeff; 

{ 

extern  double  var2[50][50],ave2[20],d; 
int  kj,i; 

double  alphau,alphal,xl[50][50],x2[50][50],xt[50][50],w; 
double  create(),standard(); 
for(k=  1  ;k<=novar;++k) 

{ 

x  1  [  1  ][k]=cur[k]+alpha*s[k]; 

} 


for(k=  1  ;k<-novar;++k) 

{ 


x2[l][k]=cur[k]; 


) 
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j— 0; 

alphau  =  alpha; 

alphal=0.0; 

while(j<=0) 


alpha=0.5*(alphau+alphal); 

for(k=l;k<=novar;++k) 

{ 

xtr[l][k]=cur[k]+s[k]*alpha-ave2[kl; 

} 

trans(xtr,xt,novar): 
mm(  b,xtr,var2, 1  ,novar,novar); 
mm(a,b,xt,  1  .novar,  1 ); 
if(a[  1  ]( 1  l>=d+.0()01 ) 

{ 

for(k=l  ;k<=novar;++k) 


alphau-alpha; 


if(a[  1  ][  1  ]<=d-.0001) 

{ 

for(k=  1  ;k<=novar;++k) 

( 

alpha! =alpha; 


/*  condition  for  exiting  the  while  loop  */ 
if(a[  1  ][  1  )<=d+.0001  &&  a(  1  ][  1  ]>=d-.0001 ) 

{ 

for(k=  1  ;k<=novar;++k) 

{ 

xtr[  1  ][k]=xtr[  1  ]  [k  J ; 

} 

j=i; 

) 

if(alphau-alphal<=. 00001 )  j=l ; 

) 

retnm(alpha); 

} 

/****»*************************************************************/ 

/*  lines  is  a  line  search  for  max  value  on  a  line  between  two  bounds  */ 

/*  this  is  done  with  secant  search  */ 

/******************************************************************/ 
double  linesd, temp, slope) 
double  temp, 1, slope; 

{ 

intj; 

double  bound  1  ,bound2,q,abs(),z,agrad(),d,e,p,rgrad,lgrad; 
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j=0; 

bound  l=temp; 

bound2=l; 

rgrad=slope; 

lgrad=agrad(l); 

if(lgrad<=0.0) 

{ 

temp  =0.0; 
goto  out; 

} 

while(j==0) 

{ 

d=rgrad; 

e=rgrad-lgrad; 

q=(temp-l); 

z=(temp-d/e*q); 

if(z  <  bound2  II  z  >  bound  1 ) 

{ 

printf("there  is  a  mistaJte\n"); 

} 

p=agrad(z); 

/*printf("this  is  p,z,l,r  %f  7c  f  7c f  %f\n”,p,z,l,temp);*/ 
if(temp  -1<=  0.00001) 

{ 

j=i; 

} 

if(p<=.000001  &&  p>=0.0) 

I 

j=i; 

) 

eise 

{ 

if(p>=0.0) 

{ 

l=z; 

lgrad=p; 

} 

else 

{ 

temp=z; 

rgrad=p; 

} 


) 

} 

temp=z; 


out:; 

retum(temp); 

} 


/** ****************************************************** ********/ 

/*  modelg  evaluates  the  model  at  x  */ 

^* ****************************************************** ***********/ 


67 


double  modelg(x) 
double  x[20|; 


{ 

extern  double  c[30],g[2()|; 
extern  int  ncoeff.novar; 
double  y,pow(),big[50][50],standard(); 
int  i; 

/♦create  evaluate  model  */ 
y=0.0; 

for(i=l;i<=novar;++i) 

{ 

big[  1  If  i]=x[i); 

} 

create(big,g,l); 

standard(big); 

for(i=l;i<=ncoeff;++i) 

{ 

y=y+c[i]*(big[l][i]); 

} 

retum(v); 

} 

I *  ***********************************************  * j 

/*  abss  determines  the  absolute  value  of  s  */ 

/*  **********************************************  *y 

double  abss(s,novar) 
double  s[20]; 
int  novar; 

{ 

double  pow(),result; 
int  i; 

result=0.0; 

for  (i=l;i<=novar;++i) 

{ 

result=pow(s[i],2.0)+result; 

) 

result=pow(result,0.5); 

retum(result); 

) 

y*  ************************************************* *********y 

/*  bk  prints  results  into  file  */ 

y***********************************************************y 

double  bk(yhat,p,tck,y) 
double  yhat.tck.y; 
int  p; 

{ 

double  up.down; 

FILE  *sp; 
up=yhat+tck; 
down=yhat-tck; 
sp=fopen("tp'7'a"); 

fprintf(sp,"  %d  %f  %f  %f  %f\n",p-l,up,yhat,down,y); 
fclose(sp); 
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} 

/*************************************************%*************:|y 

/*  this  subroutine  does  weighted  regression  */ 

^****4^***************************************%******%***********^ 

double  reg(nodat,novar,xs,var.sigma,xr,c) 
int  nodat.novar, 

double  var[50][50],xr[20],c[30],xs[10],*sigma; 

{ 

extern  int  ncoeff,flag9; 

extern  double  ave[20],sig2[20],ave2[20],var2[50][50]; 

double  mm(),y[50][50],v[50][50],xn[50],fr,pt[5],x[50][50],xt[50][50]; 

double  fabs(),ipvt[50],det[3],work[50],a[50][50],b[50][50],g; 

double  exp(),pow(),proj[50][50],avel[20],proj2[50][50]; 

int  k,l,info,counter,j,i,junk; 

double  dgefa_(),dgedi_(),create(),z[50][50],s[50][50]; 

FILE  *sp; 


lit***************************************************************  */ 

/*  beginning  of  program  to  do  local  regression  */ 

/*  *i 

j+ ********************************************************************y 


/*  read  in  data  */ 
sp=fopen("exper.dat","r"); 
printf("%d  %d\n",nodat,novar); 
t*  create  x  and  xt  matrix  */ 
for(  i=  1 ;  i<=nodat- 1 ;  ++i) 

{ 

fscanf(sp,"%f',&y[i][l]); 

for(j=l;j<=novar;++j) 

{ 

fscanf(sp,"%f',&x[i][j]); 

} 

} 

fscanf(sp,"%f',&fr); 
for(j= 1  ;j<=novar,++j) 

{ 

fscanf(sp,"%f',&xr[j]); 

} 

fclose(sp); 

y[nodatj[l]=fr, 

for(J = 1  ;j  <=novar,  ++j ) 

{ 

x[nodat][j]=xr[j]; 

} 

/*  for  each  data  pt  find  weight  and  its  contribution  to  sum  of  squares  */ 
/*  based  on  normalized  distance  from  reference  point  */ 


for<i=  1  ;i<=nodat- 1  ;++i) 

{ 

for(j=  1  ;j<=nodat;++j) 

{ 

v[i]U)=0-0; 

} 
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} 

for(i=  1  ;i<=nodat;++i) 

{ 

xn[i]=0.0; 

for(j=  1  ;j<=novar,++j) 

{ 

^xn[i]=xn[i]+pow((x[i][j]-xr[j])/xs[j],2.0); 

v[i][i]=exp(-xn[i]/2.0/novar): 

} 

/***********************************************************************/ 
/*  this  pan  is  for  distance  constraint  */ 

^* ********************************************************************** ^ 

/*  S  =  {X[I  -  l/n(ll')]X)/(n-  1)  */ 

/*  where  lis  a  vector  of  ones  */ 

for(i=  1  ;i<=nodat;++i) 

{ 

for(j = 1  ;j<=nodat; ++j) 

{ 

proj2[i]  □]=0.0- 1 .0/(nodat); 
if(i==j)  proj2[i][j]=l .0- 1 .0/(nodat); 

} 

} 

for(i=l;i<=novar,++i) 

{ 

for(j=l;j<=nodat;++j) 

{ 

xt[i][j]=x[j][i]; 

} 

} 

mm(a,xt,proj2,novar,nodat,nodat); 
mm(var2,a,x,novar,nodat,novar); 
for(j = 1  ;j  <=novar, ++j ) 

{ 

for(i=  1  ;i<=novar,++i) 

{ 

var2[i][j]=var2[i]  [j]/(nodat- 1 .0); 

} 

} 

for(i= 1  ;i<=novar,++i) 

{ 

for(j=  1  ;j<=novar,++j) 

{ 

var2[i-l][j-l]=var2[i][j]; 

} 

} 

junk=novar; 

dgefa_(var2,&junk,&junk,ipvt,&info); 
dgedi_(var2,&junk)&junk,ipvt,det,work,l  1 ); 
for(i=ncoeff;i>=0;--i) 

{ 

for(j=ncoeff;j>=0;--j) 

( 
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var2[i+l][j+l]=var2[i][j]; 

} 

/*  find  mean  vector  of  X  */ 

for(i=l;i<=novar;++i) 

{ 

ave2[i]=0.0; 

} 

for(j= 1  ;j<=nodat;++j) 

{ 

for(i=l;i<=novar,++i) 

{ 

ave2[i]=ave2[i]+x[j][i]; 

} 

} 

for(i=  1  ;i<=novar,++i) 

{ 

ave2[i]=ave2[i]/(nodat); 

} 

for(i=l;i<=nodat;++i) 

{ 

create(x,xr,i); 

} 

/*  create  projection  matrix  =  identity  matrix  */ 
for(i=l;i<=nodat;++i) 

{ 

for  (j = 1 ;  j  <=nodat; ++j ) 

{ 

proj[i][j]=0.0; 
if(i==j)  proj[i][j]=1.0; 

} 

) 

/*  find  mean  vector  of  predictor  variables  */ 

for(i=l  ;i<=ncoeff;++i) 


ave[i]=0.0; 

} 

for(j=  1  ;j<=nodat;  ++j ) 

{ 

for(i=l  ;i<=ncoeff;++i) 

{ 

ave[i]=ave[i]+x[j][i]; 

} 

} 

for(i=  1 ;  i  <=ncoeff;  ++i) 

{ 

ave[i]=ave[i]/(nodat); 

} 

mm(z,proj,x,nodat,nodat,ncoeff); 


*  4c  **  *  *  *****  *  *  ***  **  *  **  **  **  **  **  *  *  *  *  *  *  **  ****  *  *  *****  **  ***  *********** 
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/*  create  the  x  transpose  (xt)  matrix  */ 

/t*****************************************************., ,***********/ 

for(i=l;i<=ncoeff;++i) 

{ 

for(j = 1  ;j  <=nodat;  ++j ) 

{ 

xt[i][j]=z[j][i]; 

} 

} 


/*  scale  predictor  variables  by  their  average  value  */ 

/*  (this  is  done  to  aid  in  the  matrix  inversion)  */ 


for(i=l;i<=nodat;++i) 

{ 

for(j=l;j<=ncoeff;++j) 

{ 

z[i]U]=z[i]Lj]/ave[j]; 

xt[j][i]=z[i][j]; 

} 

} 

/**  ***************************************]|>************************y 

/*  b=(xt*v*x)**-l*(xt*v*y)  */ 

/*>i^**************>i^***********************************************/ 

mm(b,xt,v,ncoeff,nodat,nodat); 

mm(a,b,z,ncoeff,nodat,ncoeff); 

/*  copy  a  */ 
for(i=l  ;i<=ncoeff;++i) 

{ 

for(j=l  ;j<=ncoeff;++j) 

( 

var[i][j]=a[i]|j]; 

} 

} 

/*  find  inverse  */ 
for(i=l;i<=ncoeff;++i) 

{ 

for(j=  1  ;j<=ncoeff;++j) 

{ 

a[i-l](j*l]=a[i)[j]; 

) 

} 

junk=ncoeff; 

dgcfa_(a,&junk,&junk,ipvt,&info); 
dgedi_(a,&junk,&junk,ipvt,det,work,  1 1 ); 
for(i=ncoeff;i>=0;--i) 

{ 

for(j=ncoeff;j>=0;--j) 

{ 

a[i+l]D+l]=a[i][j]; 

} 

} 
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mm(b,a,var,ncoeff,ncoeff,ncoeff); 

y*  *************************************************************  * J 

/*  calculate  X'V**-1X  for  use  in  prediction  error  calculations  *1 

y*  ************************************************************** *y 

for(i=  1  ;i<=ncoeff;++i) 

{ 

for(j=  1  ;j<=ncoeff;++j) 

l 

var[i][j]=a[i][j]; 

} 

} 

mm(b,a,xt,ncoeff,ncoeff,nodat); 
mm(a,b,v,ncoeff,nodat,nodat); 
mm(b,a,y,ncoeff,nodat,  1 ); 
printfC'the  coefficients  are  : "); 
c[0]=fr, 

for(i=  1  ;i<=ncoeff;++i) 

{ 

printff  %f  \n",b[i][l]); 
c[i]=b[i][l]; 

printff'beta  revised  =  %f\n",c[i]/ave[i]); 

} 

printf("\n"); 

y*  *********************************************************  *  j 

/*  find  weighted  variance  */ 

y*  **************************************************** *******y 

/*  find  est  of  y  */ 
mm(a,z,b,nodat,ncoeff,  1 ); 

/*  sp=fopen("resid","w");*/ 
for(k=  1  ;k<=nodat;++k) 

{ 

z[k][l]=y[k][l]-a[k][l]; 

/*  printf("%f\n",z[k][l]); 

fprintf(sp,”%d  %f\n",k,z[k][  1  ]);*/ 

xt[l][k]=z[k][l]; 

) 

fclose(sp); 

mm(a,xt,v,  1  ,nodat,nodat); 

mm(b,a,z,  1  .nodat,  1 ); 

b[  1  ]  [  1  ] =b[  1  ]  [  1  ]/(nodat-ncoeff); 

*sigma=b[l][l]; 
printf("sigma=  %f\n",b[l][l]); 

/*  for(i=l;i<=ncoeff;-H-i) 

{ 

for(j=l  ;j<=ncoeff;++j) 

{ 

var[i][j]=b[  1  ][  1  ]*var[i][j]; 

} 

} 

*1 

} 
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/ft******************************************************/ 

/*  nran  generates  normal  random  variables  */ 

j*  **  **  *****  ***  *  ***********  *  *  *****  ****  *  ***m  *  **  *  **  *  *  4c  *  *  ***  j 

double  nran() 

{ 

int  srand().rand(); 
double  zl,rl,r2,z2; 


rl=  rand()/(pow(2.0,31.)-l): 
r2=rand()/(pow(2.,3 1 .)- 1  •); 

/*  printf("rl=%f  r2=  %f  ”,rl.r2);*/ 
zl=pow(-2.0*log(rl),0.5)*cos(2.0*3.14159*r2); 
z2=pow(-2.0*log(rl),0.5)*sin(2.0*3.14159*r2); 

/*  printf("%f  %f\n",zl,z2 );*/ 
return  (zl); 

}* 

^* ************************************************/ 
/*  vran  can  superimpose  a  random  variable  */ 

/*  on  each  of  the  control  variables  */ 

I************************************************/ 

double  vran  (cur,  curm) 
double  cur[20],curm[20]; 


{ 

double  nran(); 
int  i; 

curm[l]=nran(); 

printfC’rand  variable  =  %f\n",curm[l]); 
for(i=l;i<=5;++i) 

{ 

curm[i]=cur[i]; 

) 

} 

I*  ****************************************************************** ^ 

/*  find  projection  on  to  the  nonlinear  constraints  */ 

^******** *********************************************************** ^ 
double  w(alpha.proj) 
double  alpha,proj[20]; 


/*  find  the  projection  of  the  point  on  to  the  nonlinear  constraint  */ 

{ 

extern  double  kl,cur[20],s[20],t,aa[50][50],d; 
extern  int  pointer,novar,ncoeff,flag,pl,flag2,flag3; 
double  f[20],a[50][50],aat[50][50],aaa[50][50],h[20],proj2[20]; 
int  i  j,k,m,info,count; 

double  tk,conck(),ipvt[50] ,work[ 50] ,det[3] ,gradh(),dck,dconck(),tck; 
for(i=0;i<=pointer;++i) 

{ 

h[i]=0.0; 

} 

for(i=  1  ;i<=novar,++i) 
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( 

proj[i]=cur[i]+alpha*s[i]; 

} 

count=0; 

if(flag3==0  &.&  flag=0)  goto  q; 
wieder:; 
j=0; 

count=count+l; 
if  (count>500) 

{ 

flag2=0; 
goto  q; 

} 

if  (flag==l) 

{ 

h[pointer-flag3]=t-conck(proj); 

} 

if  (flag3==l) 

{ 

h[pointer]  =d-dconck(proj ) ; 

} 

for(i=poimer-flag-flag3+l;i<=pointer,++i) 

{ 

if(h[i]<=0.0-0.001) 

{ 

j=i; 

} 


} 

if(j==0)  goto  q; 

/*  case  where  the  first  constraint  is  binding  */ 

if(flag==l) 

{ 


if  (flag2==0)  goto  q; 
j=0; 

/*  update  projecting  matrix  */ 

gradh(f,proj); 

for(i=  1  ;i<=novar,++i) 

{ 

aa[pointer-flag3][i]=0.0-fTi]; 

} 


/*  case  where  the  second  constraint  is  binding  */ 

if(flag3=l) 

{ 

h  [pointer]  =d-dconck(proj ) ; 
if  (flag2==0)  goto  q; 

j=0; 

f*  update  projecting  matrix  */ 
dgradh(f.proj); 
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for(i=  1  ;i<=novar;++i) 

{ 

aa[pointer][i]=0.0-f{i]; 

} 

} 

for(k=  1  ;k<=pointer;++k) 

{ 

for(i=l;i<=novar,++i) 

{ 

aat[i][k]=aa[k][i]; 

} 

} 

mm(a,aa,aat,pointer,novar, pointer); 

/*  invert  matrix  */ 
for(i=  1  ;i<=pointer;++i) 

{ 

for(m=l  ;m<=pointer;++m) 

{ 

a[i-l][m-l]=a[i][m]; 

} 

} 

if  (a[0][0]==0.0) 

{ 

flag2=0; 
goto  q; 

} 

dgefa_(a,&pointer,&pointer,ipvt,&info); 
dgedi_(a,&pointer,&pointer,ipvt,det,work,01); 
for(i=pointer- 1  ;i>=0;--i) 

{ 

for(m=pointer- 1  ;m>=0;--m) 

{ 

a[i+l][m+l]=a[i][m]; 

} 

} 

mm(aaa,aat,a,novar, pointer, pointer); 

mvm(proj2,aaa,h,novar,  pointer); 
for(i=l  ;i<=novar,++i) 

{ 

proj[i]=proj[i]-proj2[t]; 

} 

goto  wieder, 
tk=conck(proj); 

/*  printf("this  is  w  end  tk  =  %fNn",tk);*/ 

q:; 


/********  ***********************  ******************* Mir************/ 

/*  wbound  returns  a  feasible  step  size  for  the  case  */ 

/*  where  at  least  one  of  the  nonlinear  constraints  is  binding  */ 

/***************************************************************/ 
double  wbound(step.proj) 
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double  step,proj[20]; 

( 

extern  double  bc[20],con_mat[50][50],s[20],cur[20],t,d; 
extern  int  novar,ncon,flag2; 
int  i,j; 

double  u[20],alphau,alphal, alpha, w(),conck(),tck,dconck(),dck; 
alpha=step; 

/*  check  to  see  if  maximum  nonlinear  step  size  results  in  a  feasible  */ 
/*  projection  on  to  the  nonlinear  constraints  */ 

alphau=bound(ncon,novar,s,cur,con_mat, be  .alpha); 
if(alphau>=step) 

{ 

alphau=step; 

/*printf("were  changing  alphauVi");*/ 

} 

alphal=0.0; 

again:; 

flag2=l; 

w(alphau.proj); 

/*  check  to  see  if  this  pt  violates  any  constraints  */ 
mvm(u,con_mat,proj  ,ncon,novar) ; 

j=i; 

for(i=  1  ;i<=ncon;++i) 

{ 

if(u[i]-bc[i]<=-.01) 

{ 

j=0; 

} 

if(d-dconck(proj)<0.0-.0l)  j==0; 
if(t-conck(proj)<0.0-.01)  j=0; 

} 


/*  if  the  point  is  feasible  then  stop  */ 


if  (j==l  &&  flag2=l) 

{ 

alpha  =alphau; 
goto  quit; 

} 

else 

/*  if  the  point  is  not  feasible,  have  the  step  size  */ 

{ 

printfO'were  in  a  new  case\n"); 
alphau  =0.5*alphau; 
goto  again; 

} 


quit:; 

retum(alpha); 
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} 

/*****************************  **************************/ 
/*  standard  divides  the  predictor  variables  */ 

/*  by  their  average  */ 

y*******************************************************^ 

double  standard(z) 
double  z[50][50]; 

{ 

extern  int  ncoeff; 

extern  double  a  -e[20],sig2[20]; 

int  i; 

for(i=l;i<=ncoeff;++i) 

{ 

z[l][i]=z[l][i]/ave[i]; 

} 

} 
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